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arc.  Certain  tidal  effects  are  also  included  in  the  adjustment.  However, 
each  diurnal  and  semidiurnal  constituent  considered  is  attributed  only  two 
adjustable  parameters:  a  global  amplitude  factor  and  a  global  phase  angle 
correction.  These  parameters  are  shown  to  be  essentially  uncontaminated 
by  the  geoidal  errors  or  by  the  systematic  orbital  errors.  Such  advantages 
would  not  exist  if  the  adjustment  were  made  in  terms  of  S.H.  tidal  coef¬ 
ficients  whose  a  priori  values  are  used  in  the  present  model  to  describe 
the  approximate  behavior  of  the  pertinent  diurnal  or  semidiurnal  constitu¬ 
ents.  The  above  two  parameters  represent  very  special  linear  combinations 
of  these  coefficients. 

The  inclusion  of  the  tidal  adjustment  in  the  overall  adjustment  of  satel¬ 
lite  altimetry  improves  the  geoidal  residuals  as  well  as  the  "observed" 
geoid  undulations  obtained  by  superimposing  the  former  on  the  adjusted 
smooth  geoid.  Such  "observed"  values  can  serve  to  produce  a  high-degree 
and  order  set  of  S.H.  potential  coefficients  via  integral  formulas  (not  via 
an  adjustment),  which  can  then  serve  in  predicting  the  desired  geophysical 
quantities.  The  geoidal  residuals  can  be  used  as  observations  in  a  sub¬ 
sequent,  ory^cond-phase,  adjustment  of  a  short-wavelength  oceanic  geoid 
in  terms  of  plint-mass  (P.M.)  magnitudes  as  parameters.  An  important  part 
of  the  de\|eloment  is  concerned  with  a  reformulation  of  the  second-phase 
adjustment  in  rarms  of  P.M.  parameters,  leading  to  a  banded  system  of 
normal  equation^.  Such  a  system  can  be  arrived  at  and  resolved  by  addres¬ 
sing  three  separate  tasks  in  what  can  be  called  a  "modified  Choleski 
algorithm":  1)  the  elimination  of  the  point  masses  from  an  observation 
equation  if  they  are  sufficiently  far  from  the  pertinent  observation  point, 
2)  the  special  arrangement  of  the  P.M.  parameters  in  the  adjustment  scheme, 
and  3)  the  resolution  of  the  resulting  system  through  an  adaptation  of  the 
well-known  Choleski  algorithm.  Whereas  previously  only  small  areas  could 
be  resolved  in  a  P.M.  approach,  the  oceanic  geoid  over  entire  ocean  basins 
can  now  be  adjusted  in  a  few  overlapping  strips  of  point  masses  by  virtue 
of  a  several -fold  reduction  in  the  run-time  and  core-space  requirements 
achieved  with  the  modified  Choleski  algorithm. 

In  a  separate  task,  a  model  for  the  deflection  rates  along  a  given  route 
is  developed  in  tensor  notations,  and  is  subsequently  specialized  for  the 
S.H.  and  P.M.  parameters.  It  is  pointed  out  that  the  rate  of  change  in  s 
when  proceeding  eastward  minus  the  rate  of  change  in  n  when  proceeding 
northward  yields  a  simple  relationship  which  can  serve  as  a  useful  mathe¬ 
matical  verification  at  either  level.  The  deflection  rates  as  well  as  the 
observational  modes  developed  and  described  in  the  previous  AFGL  reports 
are  further  presented  in  a  different  P.M.  model.  In  particular,  they  are 
transformed  from  the  context  of  single  point  masses  to  the  context  of  twin 
point  masses  (the  former  implies  a  single  P.M.  layer  while  the  latter  im¬ 
plies  a  double  P.M.  layer).  An  analysis  of  all  the  observational  modes 
considered  confirms  the  plausible  property  that  if  a  single  point  mass  is 
made  into  a  twin  point  mass,  the  effect  exercised  on  a  modeled  value  de¬ 
creases  much  more  rapidly  with  distance. 


_ Unclassified _ 

security  classification  of  this  pace(tf»>»a  D«i«  &ii«»« 


TABLE  OF  CONTENTS 


CHAPTER  SECTION  DESCRIPTION 


PAGE  NO. 


ABSTRACT  i 

1  INTRODUCTION  1 

2  ROLE  OF  THE  TIDAL  ADJUSTMENT  IN  COMPUTING 

THE  "OBSERVED"  GEOID  AS  THE  BASIS  FOR  A 
DETAILED  GRAVITY  FIELD  REPRESENTATION  5 

2.1  General  Discussion  5 

2.2  Effect  of  Tidal  Modeling  on  the  SEASAT 

Altimeter  Adjustmjnt  13 

2.3  Summary  and  Assessment  of  the  "Observed" 

Geoid  19 

3  LARGE  AREA  ADJUSTMENT  USING  THE  POINT 

MASS  PARAMETERS  22 

3.1  Notes  on  the  Strategy  in  Choosing 

the  Depth-Side  Ratio  of  Point  Masses  22 

3.2  Point-Mass  Algorithm  for  a  Banded  System 

of  Normal  Equation?  28 

4  MODEL  FOR  DEFLECTION  RATES  IN  TERMS  OF 

SPHERICAL  HARMONICS  AND  POINT  MASSES  44 

4.1  Mathematical  Background  45 

4 . 2  Spherical -Harmonic  Approach  50 

4.3  Approach  Using  a  Single  Layer  of 

Point  Masses  53 

5  OBSERVATION  MODES  IN  TERMS  OF  A 

DOUBLE  LAYER  OF  POINT  MASSES  61 

5.1  Adaptation  of  Five  Observational  Modes 

to  Twin  Point  Masse?  63 


terms  of  Twin  Point  Masses' 


73 


TABLE  OF  CONTENTS  (Continued) 


CHAPTER  SECTION  DESCRIPTION  PAGE  NO. 

6  CONCLUSIONS  78 

APPENDIX  VERIFICATION  OF  THE  POINT-MASS  MODEL 

FOR  THE  DEFLECTION  RATES  83 

ACKNOWLEDGEMENT  85 


REFERENCES 


86 


LIST  OF  FIGURES 


DESCRIPTION 


Illustration  of  the  "raw"  geoid 

Illustration  of  the  ‘‘observed"  geoid 
together  with  the  improvements  due  to 
orbital  and  tidal  adjustment 


1,  INTRODUCTION 


The  short-arc  model  used  at  AFGL,  the  goals  of  which  have  included 
the  computation  of  a  relatively  smooth  surface  approximating  the  geoid  over  the 
oceanic  areas,  is  concerned  with  adjusting  global  satellite  altimeter  data  in 
terms  of  a  truncated  set  of  spherical -harmonic  (S.H.)  potential  coefficients. 
The  altimeter  data  have  been  considered  as  the  main  source  of  observed 
quantities;  gravity  anomalies  and  other  sources  of  geopotential  information 
have  been  incorporated  via  the  weighted  S.H.  coefficients.  Six  weighted  state 
vector  parameters  per  orbital  arc  are  also  included  in  the  simultaneous  least- 
squares  process.  One  of  the  products  of  this  adjustment  is  a  revised  set  of 
S.H.  coefficients  which  can  be  used  to  predict  geoid  undulations,  gravity 
anomalies  and  other  quantities  related  to  the  disturbing  potential  (e.g.,  de¬ 
flections  of  the  vertical)  on  a  global  scale.  The  main  role  of  this  adjust¬ 
ment  is  to  provide  a  basis  for  a  more  detailed  representation  of  the  earth's 
gravity  field  to  be  carried  out  by  means  of  subsequent  adjustments. 

The  increasingly  high  quality  of  satellite  altimeter  measurements, 
in  conjunction  with  the  improvements  in  satellite  ephemeris,  calls  for  an  ad¬ 
justment  model  where  at  least  some  of  the  temporal  variations  in  the  ocean 
surface  should  be  taken  into  account.  The  most  important  variations  are  those 
caused  by  the  tide-generating  forces  of  the  moon  and  the  sun.  Of  these, 
eleven  long-period,  diurnal,  and  semidiurnal  tidal  effects  have  been  subjected 
to  adjustment  within  the  overall  altimetric  adjustment  as  described  in  [Blaha, 
1982].  The  initial  (approximate)  surface  tide  for  the  long-period  constituents 
has  been  expressed  through  a  relatively  simple  model  based  on  the  equilibrium 
tide  with  the  earth  deformation  included.  And  the  initial  surface  tide  for 


the  diurnal  and  semidiurnal  constituents  has  been  based  on  a  more  realistic 
tidal  model  which  takes  into  account  the  forces  of  friction  and  viscosity, 
the  self-gravitation,  the  ocean  loading  effects,  etc.,  and  where  the  ocean 
tide  has  been  expressed  with  the  aid  of  spherical -harmonic  tidal  coeffici¬ 
ents  developed  by  Estes  [1980].  The  adjustable  parameters  modifying  the  initial 
tidal  effects  have  been  the  amplitude  corrections  (for  all  the  constituents) 
and  the  phase  angle  corrections  (for  the  diurnal  and  semidiurnal  constituents). 

The  global  adjustment  results  in  a  "trend  surface"  approximating 
the  geoid  through  the  adjusted  S.H.  potential  coefficients,  and  in  the  resid¬ 
uals  representing  the  suppressed  geoidal  detail.  Due  to  the  relatively  small 
amplitude  of  the  geoidal  residuals,  the  introduction  of  the  spherical  approxi¬ 
mation  in  a  subsequent  adjustment  is  inconsequential.  In  the  previous  AFGL 
reports,  e.g.  [Blaha,  1979,  1980]  ,  one  such  adjustment  has  been  developed 
in  terms  of  point  masses.  In  this  approach,  the  location  of  the  point  masses 
is  stipulated  beforehand  and  only  the  point-mass  (P.M.)  magnitudes  act  as  ad¬ 
justable  parameters.  The  role  of  these  parameters  is  not  merely  to  accommo¬ 
date,  in  the  least-squares  sense,  the  residuals  from  the  previous  adjustment, 
but  also  to  provide  high-resolution  predictions  of  the  desired  geophysical 
quantities. 

The  inclusion  of  the  tidal  parameters  in  the  global  adjustment 
results  in  an  improved  trend  surface,  as  well  as  in  improved  residuals  which 
otherwise  would  absorb  the  tidal  variations.  This,  in  turn,  benefits  the 
P.M.  adjustment  discussed  above.  The  superimposition  of  the  geoidal  residuals 
on  the  adjusted  trend  surface  results  in  new  quantities,  "observed"  geoid 
undulations,  which  describe  the  surface  called  here  the  "observed"  geoid. 
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The  role  of  the  tidal  adjustment  in  improving  the  geoidal  residuals  and  the 
corresponding  "observed"  geoid  undulations  will  be  discussed  in  Chapter  2 
which  will  also  outline  the  importance  of  these  quantities  in  various 
geophysical  tasks. 

In  the  past,  the  P.M.  parameters  have  been  used  in  refining  the 
geoidal  surface  —  and  the  knowledge  of  the  earth's  gravity  field  in  general  - 
only  over  limited  areas.  This  has  been  Imputable  to  computer  limitations, 
which  are  already  felt  when  the  number  of  P.M.  parameters  surpasses  200.  In 
such  adjustments  all  the  point  masses  have  been  considered  in  conjunction 
with  every  observation  point,  i.e.,  all  the  P.M.  parameters  have  been  included 
in  every  observation  equation.  Although  this  is  the  proper  approach  from  the 
theoretical  standpoint,  it  makes  a  large-area  or  a  global  adjustment  impracti¬ 
cal  in  that  exceedingly  many  small-area  adjustments  would  have  to  be  carried 
out  in  adjacent  blocks  in  order  to  produce,  for  example,  a  geoidal  map  of  a 
large  ocean  basin.  Nonetheless,  a  simple  elimination  of  point  masses  rela¬ 
tively  far  from  an  observation  point  could  not  alone  alleviate  this  problem. 
The  solution  to  this  problem  will  be  developed  in  Chapter  3  featuring  a  string 
of  procedures  assembled  in  what  may  be  called  a  "modified  Choleski  algorithm". 

Chapter  4,  together  with  the  Appendix,  will  be  devoted  to  developing 
a  new  observational  mode  in  terms  of  both  the  S.H.  and  P.M.  parameters.  This 
mode  represents  the  rates  in  the  deflections  of  the  vertical  along  a  given 
route.  The  adjustment  model  for  the  deflection  rates  has  been  inspired  by 
the  development  of  inertial  instrumentation  used  with  growing  success  in 
various  geodetic  tasks,  notab  •  ’r»  In  .  ial  navigation. 


The  adjustment  models  in  terms  of  the  P.M.  parameters  described 
in  previous  AF6L  reports  have  been  concerned  with  a  single  P.M.  layer,  i.e., 
with  the  point  masses  all  located  at  the  same  depth.  In  Chapter  5  a  concept 
of  "twin”  point  masses  located  along  the  same  normal  to  the  earth's  surface 
will  be  developed,  corresponding  to  the  notion  of  a  double  P.M.  layer.  All 
the  observational  modes  treated  previously  will  be  given  a  twin  P.M.  formu¬ 
lation,  including  the  deflection  rates  introduced  above. 
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.  ROLE  OF  THE  TIDAL  ADJUSTMENT  IN  COMPUTING  THE  "OBSERVED" GEOI D 
AS  THE  BASIS  FOR  A  DETAILED  GRAVITY  FIELD  REPRESENTATION 


2.1  General  Discussion 

The  parameters  involved  in  the  short-arc  adjustment  of  satellite 
altimetry  are  divided  into  three  groups:  1)  corrections  to  the  spherical- 
harmonic  (S.H. )  potential  coefficients,  2)  tidal  parameters,  and  3)  cor¬ 
rections  to  the  state  vector  (s.v.)  components.  The  first  group  consists 
of  225  parameters  in  the  case  a  (14,14)  truncated  S.H.  model  is  used.  The 
second  group,  absent  in  previous  adjustments  of  satellite  altimetry  at  AFGL, 
comprises  two  parameters  per  tidal  constituent  except  for  the  long-period 
effects  where  this  number  is  one  per  constituent.  The  first  two  groups  are 
assigned  permanent  storage  in  the  computer  core.  The  third  group  consists  of 
six  s.v.  parameters  per  (short)  orbital  arc.  These  parameters  as  well  as 
the  appropriate  portions  of  normal  equations  are  assigned  reusable  storage 
which  is  one  of  the  main  features  of  the  short-arc  algorithm  described  in 
detail  in  several  AFGL  reports  and  papers  (for  example,  its  brief  review  in¬ 
cluding  several  recent  features  can  be  found  in  [Blaha,  1981]).  With  the 
aid  of  this  algorithn,  the  s.v.  parameters  are  eliminated  from  the  normal 
equations  and  are  solved  for  later,  after  the  solution  of  the  first  two 
groups  of  parameters,  arc  by  arc. 

At  the  outset  it  is  useful  to  briefly  recapitulate  the  highlights 
of  the  tidal  adjustment  as  conceived  in  [Blaha,  1982],  with  emphasis  on  diurnal 
and  semidiurnal  constituents.  It  is  based  on  a  priori  tidal  information  sup¬ 
plied  by  means  of  two  (m,m)  sets  of  S.H.  tidal  coefficients  per  constituent. 


or  2(m+l)2  coefficients.  However,  these  coefficients  are  not  being  adjusted 
at  all.  Only  two  quantities,  which  are  two  specific  functions  of  these 
coefficients,  are  adjustable  for  the  chosen  constituents:  a  parameter 
representing  a  proportional  change  in  amplitude  for  the  whole  globe,  and  a 
parameter  representing  a  change  in  phase  angle  for  the  whole  globe.  Accord¬ 
ingly,  these  parameters  do  not  change  with  location  or  time.  They  also  have 
a  plausible  physical  meaning.  In  considering  a  global  co-range  map,  the  first 
would  amount  to  changing  the  co-range  number  by  the  same  proportion  every¬ 
where;  and  in  considering  a  global  co-tidal  map,  the  second  would  amount  to 
changing  the  phase  angles  by  the  same  amount  everywhere. 

The  above  two  parameters  are  not  completely  free  to  adjust  but 
should  be  weighted  according  to  their  reliability.  Clearly,  since  the  tidal 
coefficients  have  been  obtained  by  utilizing  independent  information  subject 
to  the  solution  of  the  Laplace  Tidal  Equations,  their  weight  should  not  be 
zero.  (These  equations  take  into  account  the  forces  of  friction  and  viscosity, 
the  Coriolis  force,  the  spatial  distribution  of  depth,  as  well  as  the  self¬ 
gravitation  and  the  ocean  loading  effects  as  described  in  more  detail  in 
[Blaha,  1982],  where  heavy  use  was  made  of  [Vam'cek,  1980],  [Estes,  1980], 
[Parke  and  H endecs ho tt,  1980]  and  [Schwiderski ,  1980].)  And  since  the  above 
two  parameters  represent  certain  linear  combinations  of  the  tidal  coef¬ 
ficients,  their  weight  should  not  be  zero,  either.  Due  to  the  lack  of 
better  information,  the  "one  sigma"  assigned  for  the  independent  weighting 
of  these  parameters  in  [Blaha,  1982]  corresponded  to  0.5  for  the  relative 
amplitude  and  10°  for  the  phase  angle. 
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However,  exceptions  to  any  such  weighting  play  an  important  role, 
depending  on  the  satellite  orbital  characteristics.  For  example,  care  should 
be  exercised  when  weighting  the  constituents  K  ,  P1,  S2  and  K2  in  conjunction 
with  SEASAT  altimeter  data.  As  stated  in  [TOPEX,  1981],  in  this  case  the 
diurnal  constituents  K  and  P  are  aliased  to  a  six-month  and  a  constant 
constituent,  respectively.  The  two  semidiurnal  constituents  are  similarly 
aliased  to  a  six-month  and  a  three-month  constituent.  Since  the  useful  life 
span  of  SEASAT  amounted  to  only  about  three  months,  most  of  these  constituents 
cannot  be  properly  resolved  (even  with  aliasing  taken  into  consideration),  as 
is  already  the  case  with  the  semiannual  constituent  SSa.  Accordingly,  the 
a  priori  sigmas  associated  with  the  eight  pertinent  parameters  (two  para¬ 
meters  for  eachof  the  K  ,  P  ,  S2  and  K2  constituents)  have  been  substantially 
lowered  in  a  recent  tidal  adjustment  at  AFGL. 

A  crucial  part  of  any  discussion  involving  a  tidal  adjustment  is 
whether  or  not  the  tidal  parameters  are  contaminated  by  geoidal  and  orbital 
errors.  In  the  present  case,  the  effect  of  the  two  tidal  parameters  on  the 
modeled  sea  surface  varies  in  time.  For  example,  at  a  certain  time  when  the 
phase  is  "zero"  the  effect  of  the  amplitude  and  thus  also  of  the  amplitude 
change,  whatever  it  may  be,  on  the  geocentric  distance  to  the  sea  surface  is 
zero.  At  another  time,  the  effect  is  maximum.  The  geoid,  by  definition,  is 
time  invariant.  This  is  true  not  only  with  respect  to  the  geoidal  features 
that  are  known  and  expressed  by  some  means  (e.g.  via  S.H.  potential  coef¬ 
ficients)  but  also  with  respect  to  the  features  not  known  or  expressed,  as 
well  as  to  the  features  expressed  erroneously.  From  the  foregoing  it  follows 
that  the  present  two-parameter  solution  should  not  be  influenced,  at  least 
in  theory,  by  the  presence  of  geoidal  errors. 
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A  different  question  arises  if  one  considers  the  observational 
noise,  assumed  random.  If  the  sigma  (standard  error)  of  this  noise  is 
relatively  large,  many  of  the  local  tidal  features  may  not  be  recoverable 
due  to  numerical  difficulties.  But  the  present  two  parameters  pertain  to 
all  of  the  features,  both  large  and  small,  in  all  of  the  oceans,  associated 
with  the  given  constituent.  Thus,  changes  in  these  parameters  affect  the 
least-squares  sum  more  than  could  be  attributed  to  the  comparable  changes  in 
the  individual  tidal  coefficients  and  their  various  combinations  (especially 
the  combinations  giving  rise  to  small  local  tidal  features).  Accordingly, 
these  parameters  can  be  resolved  better  than  other  tidal  quantities  could. 
However,  a  good  quality  in  the  solution  can  be  achieved  only  if  these  tidal 
parameters  are  sufficiently  independent  of  the  systematic  errors  in  the  ad¬ 
justment.  The  independence  with  regard  to  geoidal  errors  has  been  shown  in 
the  preceding  paragraph. 

The  effect  of  orbital  errors  on  the  recoverability  of  tidal  para¬ 
meters  causes,  in  general,  the  greatest  concern.  In  this  context,  the 
reference  heavily  relied  upon  will  be  [Estes,  1980],  Its  opening  statements 
acknowledge  that  the  efforts  to  extract  ocean  tide  information  from  satellite 
altimetry  data  have  been  largely  unsuccessful,  mostly  due  to  errors  in  the 
orbit  determination.  In  particular,  it  states: 

"The  tidal  frequencies  are  well  determined  from  astronomical 
considerations,  and  computer  simulations  by  Zetler  and  Maul 
(1971)  and  by  Won  et  al .  (1977)  strongly  indicate  the  pos¬ 
sibility  of  recovering  the  principal  tides  with  widely  spaced 
time  sampling  using  this  type  of  data.  However,  these  simula¬ 
tions  have  treated  measurement  noise  and  random  orbit  error, 
but  have  not  considered  the  aliasing  effects  of  systematic 
orbit  error,  which  over  local  regions  can  resemble  a  tide  signal." 
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At  this  point  one  notices  the  term  "local  regions",  which  can  be 
contrasted  to  the  global  character  of  the  above  two  parameters.  The  same 
reference  further  states  that  unmodeled  systematic  errors  are  significant 
and  must  be  considered  in  the  analysis  for  tidal  structure  from  satellite 
altimeter  data.  Here  the  term  "tidal  structure"  clearly  applies  to  all  of 
the  features  pertaining  to  a  given  tidal  constituent,  not  merely  to  the 
amplitude  and  phase  changes  made  uniform  for  the  whole  globe. 

The  main  result  of  the  above  reference  is  expressed  on  page  92: 

"The  present  results  demonstrate  that  with  adequate  modeling 
of  the  ephemeris  errors  to  reduce  the  aliasing  of  the  tides, 
simultaneous  estimation  of  the  tide  model  and  radial  orbit 
error  model  parameters  will  accurately  recover  the  principal 
tidal  constituents." 

It  will  be  shown  that  the  present  model  is  different  from  Estes'  (and  from 
that  of  the  other  investigators  as  quoted  above)  and  that  it  does  not  require 
model  parameters  for  systematic  orbital  errors.  The  key  consideration  lies 
in  its  scope,  in  that  it  does  not  seek  to  recover  the  entire  tidal  consti- 
tuent  as  described  by  2(m+l)  S.H.  tidal  coefficients.  Instead,  it  attempts 
to  recover  only  two  parameters  which,  furthermore,  are  weighted. 

The  concern  addressed  in  [Estes,  1980]  is  the  separability  of  the 
orbit  error  and  the  tide  (see  p.76).  But  the  term  "tide"  implies  all  of  the 
S.H.  tidal  coefficients.  This  is  clear,  for  example,  from  p.83  where  the 
total  number  of  adjustment  parameters  is  given  as 

220xNC  +  50x(2NP+IB+IS). 

The  first  term  pertains  to  all  the  S.H.  tidal  coefficients  for  all  the  tidal 


constituents.  The  second  term  then  includes  the  model  parameters  for 
systematic  orbital  errors  on  all  of  the  50  satellite  arcs  considered.  One 
further  notices  the  statement  on  page  81: 

"In  the  least-squares  recovery  for  the  model  parameters, 
the  a  priori  values  of  the  expansion  coefficients  are  set 
to  zero." 

At  this  point,  the  S.H.  coefficients  in  [Estes,  1980]  are  seen  as: 

a)  not  known  a  priori, 

b)  given  zero  values  a  priori, 

c)  given  zero  weight  a  priori, 

d)  all  subject  to  adjustment. 

On  the  other  hand,  the  S.H.  tidal  coefficients  in  the  present  adjustment  are 

a)  known  a  priori, 

b)  given  their  a  priori  values  obtained  from  an  independent 
source, 

c)  not  considered  for  weighting,  etc.,  because  they  are  not 
adjustable  parameters  individually;  in  particular, 

d)  only  two  special  combinations  are  subject  to  adjustment, 
they  are  assigned  two  parameters  and  these  are  attributed 
weights  according  to  their  reliability. 

Except  for  the  two  special  kinds  of  changes  as  embodied  by  the 
above  two  parameters,  the  whole  configuration  (in  space  and  time)  of  the 
tidal  constituent  in  question  is  adopted  completely  from  the  a  priori 
information.  Clearly,  if  one  wanted  to  resolve  all  of  the  S.H.  tidal 


coefficients,  one  would  be  faced  with  the  necessity  of  modeling  the  orbital 
errors  as  demonstrated  by  Estes  [1980].  In  circumventing  such  a  task,  the 
present  model  represents  a  limited  adjustment  which  can  still  Improve  the 
tidal  knowledge,  reduce  the  residuals  and,  perhaps,  indicate  deficiencies  in 
the  a  priori  tidal  information  by  using  satellite  altimetry  as  an  indepen¬ 
dent  source.  The  term  "limited  adjustment"  should  be  understood  a;  as  in¬ 
corporating  only  two  parameters  per  constituent,  and  b)  as  considering  even 
these  parameters  weighted.  Due  to  the  weighting,  unless  there  exists  a 
strong  inconsistency  between  the  altimetry  and  the  a  priori  tidal  information, 
the  corrections  to  these  parameters  are  expected  to  be  relatively  small. 

The  characteristics  of  the  present  model  are  further  elucidated  if 
one  recalls  that  a  tidal  constituent  is  described  by  several  properties, 
such  as  the  number  and  locations  of  amphidromic  points,  the  shape  of  the 
co-phase  lines  (emanating  from  the  above  points)  anywhere  on  the  globe,  the 
shape  of  the  co-range  lines  anywhere  on  the  globe,  etc.  These  properties 
are  described  by  the  a  priori  information  considered  here  to  be  the  S.H. 
tidal  coefficients.  If  these  coefficients  were  subject  to  adjustment, 
changes  in  any  or  all  of  them  would  bring  changes  to  any  or  all  of  the  above 
properties  and  vice  versa.  Thus,  if  altimeter  data  were  contaminated  by 
systematic  errors  over  a  given  area,  the  coefficients  resolved  in  an  adjust¬ 
ment  would  in  general  all  be  affected  to  a  certain  extent.  This,  in  turn, 
would  affect  any  or  all  of  the  tidal  properties  in  that  area. 

However,  the  present  tidal  model  allows  for  only  two  very  special 
changes  in  two  specific  properties:  1)  the  change  in  co-range  numbers  by 
the  same  proportion  over  the  whole  globe  and  2)  the  change  in  co-phase 


-11- 


numbers  by  the  same  amount  over  the  whole  globe,  out  of  an  infinite  number 
of  all  possible  global  and  local  changes.  As  it  has  just  been  indicated, 
the  local  properties  respond  to  systematic  orbital  changes  in  various  areas 
and  during  various  time  periods.  But  only  some  very  special  systematic 
errors  in  the  orbit,  present  all  over  the  globe  all  of  the  time,  could 
introduce  changes  in  the  above  two  special  properties.  Even  though  syste¬ 
matic  orbital  errors  exist  they  tend  to  "average  out"  in  space  and  time, 
especially  periodic  errors.  One  could  conclude  these  heuristic  arguments 
by  stating  that  if  any  properties  of  a  tidal  constituent  are  insensitive 
to  such  systematic  orbital  errors,  it  is  those  representing  the  very  special 
global  shifts  in  the  co-range  and  co-phase  lines  described  by  the  two  tidal 
parameters  in  the  present  adjustment  model. 


During  recent  data  reductions  at  AFGL,  some  6,700  arcs  of  SEASAT 


altimetry  have  been  adjusted  in  the  short-arc  mode,  where  the  arcs'  durations 
have  been  seven  minutes  or  less.  This  adjustment  has  yielded  for  M2  a  22% 
reduction  in  amplitude  and  a  0.7°  change  in  phase  angle  as  given  by  the 
initial  set  of  (12,12)  S.H.  tidal  coefficients  supplied  by  Estes  [1980]. 

Such  small  changes  indicate  a  realistic  model  for  this  constituent  whose 
effect  is  twice  to  several  times  greater  than  the  effect  of  the  other  tidal 
constituents.  The  permanent  tide  is  an  exception  to  this  statement,  but  its 
effect  is  still  lower  than  that  of  M2;  mathematically,  it  can  be  related  to 
the  earth's  flattening. 

The  inclusion  of  the  tidal  effects  into  the  adjustment  has  re¬ 
sulted  in  slight  overall  reductions  in  the  size  of  the  residuals  and  the 
related  quantities,  as  compared  to  the  original  adjustment  where  these 
effects  were  not  considered.  But  to  better  understand  the  relative  small¬ 
ness  of  these  reductions  the  difference  between  the  "tide"  (denoted  by  the 
index  "t")  and  the  "no  tide"  (denoted  by  the  index  "nt")  results  are  pre¬ 
sented,  where  v  is  the  residual  and  Au  is  the  "up"  correction  to  the  s.v. 
components : 

--  ,t  -  vnt  . 

A(AU)  =  Au,  -  AU  . 

t  nt 

From  randomly  distributed  satellite  passes  all  over  the  globe,  it  has  been 
found  that  the  average  Av  and  average  A(Au)  are  nearly  zero;  their  RMS  values 
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have  been  computed  as 


RMS  (Av)  =  0.047m  , 

RMS[  A(Au)]  =  0.187m  . 

Using  a  majority  of  the  adjusted  arcs,  vt  and  vnt  have  been  found 
to  have  a  nearly  zero  average  value;  further, 

RMS  ( v  )  =  1.733m  ,  RMS  (v  J  =  1.734m  . 

t  nt 

Clearly,  the  RMS  improvement  in  v  as  compared  to  values  vnfc  is  very  small 
because  RMS  (Av),  although  itself  reaching  5cm,  is  small  by  comparison.  In 
fact,  it  turns  out  that  1.7332+  0.0472  =  1.7342,  which  could  be  given  a  rough 
interpretation  where  Av  could  represent  random  errors  attached  to  v  in  order 
to  yield  v  in  a  stochastic  process  on  a  sphere.  A  reduction  in  RMS  (Au)  is 
also  very  small,  represented  by 

RMS  (Au  )  =  1.869m  ,  RMS  (Au  J  =  1.875m  . 
t  nt 

In  analogy  to  the  above,  RMS  [a( Au)]  is  small  when  compared  to  these  values 
and  a(au)  can  be  given  an  interpretation  similar  to  that  for  Av.  The  RMS 
value  of  1.869m  (and  also  1.875m)  in  the  "up"  correction  to  the  state  vectors 
compares  very  well  with  the  a  priori  sigma  of  1.6m  associated  with  the  NSWC 
precise  ephemeris  used. 

One  of  the  most  important  results  of  SEASAT  altimetry  and  tidal 
adjustment  are  the  "observed"  geoid  undulations,  denoted  N„0„,  describing 
the  "observed"  geoid.  The  quotes  are  used  to  indicate  that  certain 
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improvements  have  taken  place  in  relation  to  more  directly  observed  undula¬ 
tions,  called  here  the  "raw"  geoid  undulations.  The  latter,  denoted  N„r„, 
describe  the  "raw"  geoid  illustrated  in  Figure  1.  A  "raw"  geoid  undulation 
is  given  essentially  as  the  initial  radial  distance  (from  the  geocenter)  to 
the  satellite,  minus  the  radial  distance  to  the  altimeter  foot-point  on 
the  ellipsoid,  minus  the  altimeter  measurement,  plus  the  smal1  correction 
"d",  due  to  the  earth's  flattening,  mentioned  e.g.  in  Section  4.1  of  [Blaha, 
1982].  The  "observed"  geoid,  depicted  in  Figure  2,  can  be  constructed  from 
the  "raw"  geoid  through  two  refinements,  1)  by  utilizing  the  adjusted 
(rather  than  initial)  s.v.  components  and  2)  by  considering  the  adjusted 
tidal  effects. 

The  direct  representation  of  the  "observed"  geoid  is  made  through 
the  first  equation  in  Figure  2, 

N„0„  =  Na  +  v  , 

where  Na  is  the  adjusted  geoid  undulation  computed  from  the  adjusted  set  of 
S.H.  potential  coefficients  and  v  is  the  residual  as  introduced  earlier. 

The  representation  displaying  the  difference  between  the  "observed"  and 
the  "raw"  geoids  is  given  in  the  second  equation  of  the  same  figure  as 

N,|g||  _  N,|r„  ”  ^2  * 

where  -1  is  the  improvement  due  to  the  orbital  adjustment  (i.e.,  adjustment 
of  the  s.v.  components)  and  -i  is  the  improvement  due  to  the  adjusted  tidal 
effects.  At  a  mid-arc  epoch  the  quantity  -1  coincides  with  the  "up"  cor¬ 
rection  whose  characteristics  have  been  already  discussed.  The  RMS  of  -1 
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N„rii  *  N°  +  L  .  N,.r„  *  "raw"  geoid  undulation 


N°  =  geoid  undulation  correspond¬ 
ing  to  approximate  spherical- 
harmonic  potential  coefficients 


L  =  constant  term  of  altimeter 
observation  equation  in  the 
spherical -harmonic  model 


Figure  1 

Illustration  of  the  "raw"  geoid 


Altimeter  observation  from 
adjusted  satellite  position 


N,,o"  =  N3  +  v  «••••  Nn()„  *  "observed"  geoid  undulation 

Na  *  geoid  undulation  corresponding  to 
adjusted  S.H,  potential  coefficients 

v  ■  residual  in  the  S.H.  model 


Nv,  =  NnrM  -  ij  -  i-2  .  -i-j  =  improvement  due  to 

orbital  adjustment 


-i?  =  improvement  due  to 
tidal  adjustment 


Figure  2 

Illustration  of  the  "observed"  geoid  together  with  the 
improvements  due  to  orbital  and  tidal  adjustment 
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has  been  computed  from  the  randomly  distributed  arcs  mentioned  earlier  as  29cm, 
which  is  indeed  a  realistic  value  for  open  ocean. 

The  improvements  in  Nllo„  due  to  the  tidal  adjustment  can  be  assessed 
by  comparing  the  corresponding  results  with  the  non-tidal  adjustment  according 
to 

ANll0„  *  AN3  +  Av  , 

where  Av  was  already  defined  and  where 


ANa 

=  Na  - 

Na  , 

t 

nt 

ANllQ.. 

-  N"o 

"t  '  N"0"nt 

From  the  same  arcs  as  before,  the  RMS  values  have  been  computed  as 

RMS  (ANa)  =  0.111m  , 

RMS  (AN,,  ,,)  -  0.121m  . 

The  improvements  in  Na  and,  more  importantly,  in  N„0„  due  to  the  tidal  adjustment 
are  thus  seen  to  be  at  the  decimeter  level. 

In  the  past,  the  "observed"  geoid  undulations  along  the  SEASAT  passes 
were  useful  in  a  study  of  ocean-bottom  phenomena,  especially  the  bathymetry. 

3ut  the  benefits  of  the  "observed"  geoid  do  not  end  there.  For  example,  after 
having  a  smooth  "trend  surface"  (corresponding  to  a  14,14  S.H.  model  or  other¬ 
wise)  removed  from  this  geoid,  the  resulting  differenced  undulations  can 
serve  in  expressing  the  geoidal  surface  to  a  desired  resolution,  whether  on  a 
regional  or  a  global  scale. 
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2.3  Summary  and  Assessment  of  the  "Observed"  Geoid 


The  adjustment  model  of  satellite  altimetry  presented  herein  has 
been  based  on  the  short-arc  algorithm.  In  addition  to  the  newly  introduced 
tidal  parameters  the  adjustment  model  includes  a  set  of  low  degree  and  order 
spherical -harmonic  potential  coefficients,  typically  a  (14,14)  set,  and  a 
set  of  orbital  parameters  comprising  six  state  vector  components  per  orbital 
arc.  All  these  parameters  are  subject  to  a  simultaneous  least-squares  ad¬ 
justment.  In  earlier  as  well  as  recent  adjustments,  the  potential  coef¬ 
ficients  and  the  state  vector  components  have  been  weighted  according  to 
their  reliability  specified  at  the  source  (GEM  10  coefficients,  NSWC  precise 
ephemeris).  The  weights  for  the  tidal  parameters  have  been  introduced  in 
[Blaha,  1982]  as  corresponding  to  the  a  priori  sigma  of  1  (long-period 
constituents)  or  0.5  (diurnal  and  semidiurnal  constituents)  for  the  relative 
amplitudes,  and  to  the  a  priori  sigma  of  10°  for  the  phase  angles.  However, 
the  weights  have  been  substantially  increased  in  the  case  of  Kx,  S2  and 
K2  constituents  which  cannot  be  properly  resolved  by  SEASAT  altimetry  and 
would  otherwise  be  aliased  into  constituents  of  lower  frequencies. 

The  adjustment  including  the  tidal  effects  represents  an  improve¬ 
ment  over  an  earlier  non-tidal  adjustment  of  altimeter  observations.  The 
differences  between  the  two  kinds  of  adjustments  translate  into  about  11cm 
RMS  for  the  adjusted  geoid  (computed  from  the  adjusted  coefficients),  5cm  RMS 
for  the  residuals,  19cm  RMS  for  the  "up"  component  of  the  state  vectors  and 
12cm  RMS  for  the  "observed"  geoid.  The  latter  is  obtained  by  adding  the 
residuals  to  the  corresponding  adjusted  geoid  undulations,  and  is  equivalent 
to  the  directly  observed  "raw"  geoid  improved  through  the  adjustment 


-19- 


corrections  to  the  state  vectors  and  through  the  adjusted  tidal  effects. 
Since  the  "observed"  geoid  plays  a  major  role  in  a  detailed  determination 
of  the  earth's  gravity  field,  an  improvement  at  a  decimeter  level  can  be 
significant. 

The  features  offered  by  the  short-arc  altimetry  adjustment  in¬ 
cluding  the  tidal  effects  are  now  briefly  summarized: 

-  The  adjustment  is  global  in  nature  and  is  capable  of  discerning 
important  tidal  effects. 

-  The  inclusion  of  the  tidal  effects  has  lead  to  a  5cm  RMS  improve¬ 
ment  in  the  residuals  (as  opposed  to  a  "non-tidal"  adjustment) 
with  the  consequence  of  a  slight  reduction  in  the  overall  magnitude 
of  the  residuals. 

-  A  similar  conclusion  has  been  reached  for  the  "up"  corrections  to 
the  state  vectors  (the  improvement  has  been  19cm  RMS),  whose  over¬ 
all  RMS  value  of  1.87m  compares  very  well  with  the  a  priori  sigma 
of  1.6m. 

-  The  improvements  in  the  adjusted  geoid  (computed  from  the  adjusted 
spherical -harmonic  potential  coefficients)  and  in  the  "observed" 
geoid  have  amounted  to  11cm  and  12cm  RMS,  respectively. 

-  The  "observed"  geoid  and  the  corresponding  residuals  can  play  an 
important  role  in  detecting  and  studying  a  variety  of  ocean-bottom 
Phenomena. 

-  The  "observed"  geoid  can  serve  in  refining  the  global  oceanic  geoid 
via  various  prediction  techniques.  One  such  technique  uses  the 
point-mass  parameters  as  explained  in  previous  AFGL  reports  as  well 


as  in  Chapters  3  and  5  herein.  As  an  example  of  another  technique, 
the  "observed"  geoid  can  be  used  to  produce  a  high  degree  and  order 
set  of  spherical -harmonic  potential  coefficients  via  integral  for¬ 
mulas,  which  can  then  serve  in  predicting  geoid  undulations  and 
other  geophysical  quantities  related  to  the  disturbing  potential 
(e.g.,  gravity  anomalies,  deflections  of  the  vertical,  etc.). 

-  One  can  conclude  that  in  taking  advantage  of  the  adjusted  state  vec¬ 
tors  and  the  adjusted  tidal  effects,  and  of  the  high-resolution 
residuals  at  the  same  time,  the  "observed"  geoid  leads  to  improvements 
in  a  detailed  representation  of  the  earth's  gravity  field.  As  one 
example  of  using  these  high-resolution  residuals  in  an  independent 
approach,  techniques  are  being  developed  at  AFGL  to  compute  gravity 
anomalies  from  the  "observed"  geoid  or  the  corresponding  residuals 
along  certain  ocean  trench  areas,  which  have  already  shown  a  good 
agreement  with  the  ground  truth. 


3,  LARGE  AREA  ADJUSTMENT  USING  THE  POINT  MASS  PARAMETERS 


3.1  Notes  on  the  Strategy  in  Choosing 
tKe  Depth-Side  Ratio  of  froint  Masses 


In  one  of  the  original  reports  on  a  point-mass  (P.M.)  adjustment, 
[Needham,  1970]  suggested  the  depth/side  ratio  for  point  masses,  d/s,  as 
0.8/1.  That  study  concentrated  to  a  great  extent  on  gravity  anomalies. 

In  the  AFGL  report  [Blaha,  1979],  both  gravity  anomalies  and  geoid  undula¬ 
tions  were  taken  into  consideration  in  the  computer  simulations  which 
eventually  lead  to  the  recommendation  that  the  above  ratio  be  doubled,  in 
the  sense  d/s=  1.6/1.  However,  due  to  the  computer  core  limitations,  both 
the  simulated  data  and  the  spherical -harmonic  (S.H.)  coefficients  recovered 
from  the  first  adjustment  corresponded  to  only  a  very  low  degree  and  order 
S.H.  expansion.  In  a  typical  case,  "errorless"  data  were  generated  with  a 
given  (10,10)  set  and  the  first  (global)  S.H.  adjustment  was  performed  in 
terms  of  a  (6.6)  set.  This  means  that  the  P.M.  parameters  were  to  accom¬ 
modate  the  residuals  from  the  first  adjustment,  which  only  contained  the 
gravity  field  information  within  (7,7)  to  (10,10)  truncation. 

At  the  present,  the  first  adjustment  alone  corresponds  to  a  (14, 
14)  set  of  S.H.  potential  coefficients  and  the  P.M.  parameters  (i.e.,  P.M. 
magnitudes)  are  distributed  typically  in  a  2°x2°  equilateral  grid  which 
serves  to  accommodate  the  residuals  within  an  up  to  about  (90,90)  equivalent 
truncated  sat.  It  is  clear  that  the  P.M.  depth  should  be  related  to  the 
information  contained  within  (15,15)  to  (90,90)  truncated  sets,  whereas  the 
simulations,  as  stated  above,  could  only  take  into  consideration  (7,7)  to 
(10,10)  sets.  For  this  reason  the  ratio  1.6/1  is  not  considered  final. 


Another  motivation  for  re-examining  the  ratio  d/s=  1.6/1  is 
practical  in  nature.  If  all  the  P.M.  parameters  were  included  in  the  for¬ 
mation  of  every  observation  equation  as  was  the  case  with  the  above  computer 
simulations,  the  computer  run-time  would  be  exceedingly  high,  especially  if 
a  detailed,  large-scale  adjustment  of  the  oceanic  geoid  were  required. 
Furthermore,  the  computer  storage  difficulties  would  be  almost  insurmountable 
without  the  use  of  an  algorithm  tailored  for  a  banded  structure  of  normal 
equations.  Such  an  algorithm,  described  in  the  next  section,  stipulates 
that  beyond  a  certain  spherical  cap  centered  on  an  observation  point  (or, 
equivalently,  beyond  a  certain  cut-off  distance)  the  P.M.  parameters  are 
ignored  in  the  pertinent  observation  equation.  This  introduces  approxima¬ 
tions  in  the  adjustment  model  which,  however,  are  inconsequential  if  the 
effect  of  the  neglected  P.M.  parameters  on  the  modeled  value  at  the  observa¬ 
tion  point  is  very  small  compared  to  the  effect  of  the  point  masses  located 
in  the  vicinity  of  this  point.  In  using  a  special  arrangement  of  the  P.M. 
parameters,  the  above  cut-off  distance  represents  a  transition  from  a  full 
matrix  of  normal  equations  (N)  to  a  banded  matrix.  The  next  paragraph 
further  elaborates  on  this  property. 

In  a  P.M.  adjustment  model  the  correlation  between  the  parameters 
decreases  as  the  distance  between  the  corresponding  point  masses  grows 
larger.  In  using  a  special  design  for  grouping  the  parameters,  most  of  the 
off-diagonal  elements  in  the  matrix  N  can  be  made  very  small  (in  the  sense 
that  the  further  from  the  main  diagonal,  the  smaller  the  elements  become). 
Without  a  significant  loss  of  accuracy  some  of  these  elements  can  be  set  to 
zero,  resulting  essentially  in  a  banded  system  of  normal  equations.  In 
practice,  this  task  can  be  accomplished  at  the  level  of  observation 
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equations,  in  conjunction  with  the  cut-off  distance  described  above.  The 
resulting  system  can  then  be  resolved  with  a  great  efficiency  by  a  modified 
Choleski  algorithm.  Such  an  algorithm  allows  for  the  inclusion  of  many 
more  P.M.  parameters  in  the  adjustment,  perhaps  ten-fold,  than  would  other¬ 
wise  be  possible  due  to  the  computer  storage  and  run-time  limitations. 

This  approach  can  be  used  in  conjunction  with  both  the  single  and  the  twin 
point  masses  (the  latter  are  described  in  Chapters).  As  a  result  of  the 
above  modeling  approximation,  a  detailed  resolution  of  the  geoidal  surface 
and  the  earth's  gravity  field  is  made  possible  on  a  large,  or  even  global, 
scale. 

The  efficiency  of  the  modified  Choleski  algorithm  increases,  both 
from  the  run-time  and  computer  core  standpoints,  if  the  bandwidth  in  the 
matrix  N  decreases.  In  considering  a  suitable  arrangement  of  the  parameters 
according  to  the  P.M.  locations,  the  bandwidth  is  linked  directly  to  the 
size  of  the  spherical  cap,  in  the  sense  that  a  smaller  cut-off  distance  re¬ 
sults  most  often  in  a  smaller  bandwidth.  But  the  cut-off  distance  can  be 
reduced  without  unduly  compromising  the  rigor  of  the  solution  only  if  the 
"kernel"  function,  representing  the  effect  of  a  point  mass  on  the  modeled 
value,  decreases  sufficiently  with  distance.  When  examining  the  relative 
values  of  this  function  at  various  n-multiples  of  the  P.M.  separation  s, 
one  could  be  guided  by  the  following  very  approximate  classifications: 

40%  -  50%  effect  remains  .  cut-off  decision  is  marginal, 

20%  -  40%  effect  remains  .  cut-off  decision  is  sufficient, 

under  20%  effect  remains  .  cut-off  decision  is  more  than 

sufficient. 
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Clearly,  a  shallower  point  mass  results  in  a  kernel  function 
which  decreases  more  rapidly  with  increasing  distance.  If  the  ratio  d/'s  = 
1.6/1  were  used,  the  cut-off  decision  would  have  to  be  made  in  favor  of 
n=  3  or  ns  4  with  regard  to  some  of  the  observational  modes  described  in 
[Blaha,  1980],  especially  the  gee  d  undulations.  Being  derived  directly 
from  satellite  altimetry,  the  geoid  undulations  represent  the  most  impor¬ 
tant  mode  in  the  P.M.  adjustment.  With  the  cut-off  distance  of  3s  or  4s 
the  computer  requirements  would  still  be  unrealistically  high  for  a  geoidal 
adjustment  in  a  large  ocean  basin  such  as  the  Indian  Ocean,  etc.  However, 
if  d/s  is  lowered  to  the  original  ratio  of  0.8/1,  the  same  approximation 
level  is  compatible  with  the  cut-off  distance  of  only  1.5s  or  2s.  One 
could  lower  this  ratio  even  further  but  it  might  not  be  desirable  for  other 
reasons  (in  [Needham,  1970],  this  ratio  represented  a  close  relationship 
between  the  P.M.  blocks  and  the  gravity  anomaly  blocks).  Be  that  as  it  may, 
the  0.8/1  ratio  adopted  in  conjunction  with  the  cut-off  distance  of  1.5s 
appears  satisfactory  for  the  purpose  at  hand,  as  is  illustrated  next. 

In  view  of  the  new  ratio  d/s=  0.8/1,  the  kernel  functions  for  the 
five  observational  modes  analysed  in  [Blaha,  1980]  must  be  re-examined.  The 
value  s  of  the  P.M.  separation  now  corresponds  to  445  km  (4°  in  arc  measure), 
whereas  the  value  of  d  remains  at  350  km.  The  theoretical  outcome  of  such 
an  analysis  is  not  tied  to  a  specific  P.M.  separation  (for  example,  2°  in 
arc  measure  could  be  used  as  well,  etc.)  because  one  is  interested  merely  in 
the  relative  P.M.  effect  as  a  function  of  n.  The  results  presented  in  the 
above  reference,  computed  for  a  point  mass  whose  magnitude  is  a  10~6th  part 
of  the  earth's  mass,  are  modified  as  follows: 
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Geo id  undulations 


n  =  1.5 


55.2  m  (47.6%) 
43.6  m  (37.6%) 
35.8  m  (30.9%) 


r* 


n  =  0 

..116  m  (100%) 

n  =  1.5 

...  55.2  m  (47.6%) 

n  =  0.5  . . . 

..  98.9  m  (85.3%) 

n  =  2 

...  43.6  m  (37.6%) 

n  =  1 

..  73.1  m  (63.0%) 

n  =  2.5  .. 

. ..  35.8  m  (30.9%) 

Gravity  anomalies 

n  =  0 

..  290  mgal  (100%) 

n  =  1 

. . .  62.0  mgal  (21.4%) 

n  =  0.5  ... 

..  172  mgal  (59.3%) 

n  =  1.5 

...  21.2  mgal  (  7.3%) 

Deflections 

of  the  vertical 

n  =  0 

...  0.0"  (  0%) 

n  =  1 

...  20.5"  (79.8%) 

n  =  0.25  .. 

...  18.0"  (70.0%) 

n  =  1.5 

...  13.2"  (51.4%) 

n  =  0.5  . . 

...  25.4"  (98.8%) 

n  «  2 

...  8.7"  (33.9%) 

n  =  0.57  .. 

...  25.7"  (100%) 

n  =  2.5 

...  6.0"  (23.3%) 

n  =  3 

...  4.4"  (17.1%) 

Horizontal 

gradients  of  gravity  disturbance 

i.  -  0 

...  0.00  E  (  0%) 

n  =  0.5 

...  7.54  E  (97.3%) 

n  =  0.25  .. 

...  6.68  E  (86.2%) 

n  =  1 

...  3.44  E  (44.4%) 

n  =  0.45  . . 

...  7.75  E  (100%) 

n  =  1.5 

...  1.33  E  (17.2%) 

n  =  2 

...  0.58  E  (  7.5%) 

Vertical  gradients  of  gravity  disturbance 


n  =  0 


n  =  1 


n  =  1.5 


n  -  0.5 


18.59  E  (100%) 
6.97  E  (37.5%) 
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0.67  E  (  3.6%) 
-0.19  E  (-1.0%) 


As  it  appears  from  the  above,  the  cut-off  distance  at  n  =  1.5  is 
more  than  sufficient  for  all  the  observational  modes  except  geoid  undula¬ 
tions  and  deflections  of  the  vertical,  for  which  it  can  be  termed  marginal. 
However,  due  to  practical  considerations  (especially  the  bandwiath)  this 
distance  is  adopted  in  all  cases.  If  the  ratio  d/s=  1.6/1  had  been  used  for 
undulations,  n=l,5  would  have  corresponded  to  the  effect  of  about  75%.  This 
would  have  lowered  the  rigor  of  the  least-squares  adjustment,  increased  the 
residuals,  etc.  A  test  adjustment  for  a  limited  area  in  the  North  Atlantic 
has  confirmed  that  with  the  same  cut-off  distance  of  1.5s,  the  0.8/1  ratio 
leads  to  a  very  slightly  lower  RMS  residual  than  the  1.6/1  ratio.  In 
particular,  the  1.6/1  ratio  has  produced  a  1.59  m  RMS  residual,  the  1.2/1 
ratio  has  produced  a  1.58  m  RMS  residual,  and  the  0.8/1  ratio  has  produced 
a  1.56  m  RMS  residual  using  the  same  data  of  SEASAT  altimetry.  The  computer 
simulations  mentioned  earlier  yielded  a  different  outcome  (a  higher  ratio 
resulted  in  a  lower  RMS  residual)  due,  to  a  great  extent,  to  the  theoretical 
treatment  in  which  the  cut-off  distance  was  infinite,  and  also  due  to  a  low 
level  of  geoidal  variation  as  implied  by  the  generated  data. 


3.2  Point-Mass  Algorithm  for  a  Banded  System  of  Normal  Equations 


The  cornerstone  of  the  present  development  is  an  application  of 
the  Choleski  algorithm  to  a  large  scale  geoid  determination  from  satellite 
altimetry.  The  reason  for  using  this  algorithm  is  that  it  lends  itself 
well  to  an  accurate  and  efficient  solution  of  nonsingular  (therefore, 
positive-definite)  systems  of  normal  equations,  especially  banded  systems. 

The  latter  are  characterized  by  a  great  number  of  zero  elements  in  the 
matrix  of  normal  equations  (N),  located  outside  of  a  central  band.  How¬ 
ever,  before  arriving  at  a  banded  system  a  reformulation  of  the  least-squares 
problem  often  has  to  take  place. 

In  the  present  context  one  strives  for  the  smallest  possible  number 
of  parameters,  i.e.,  point-mass  (P.M.)  magnitudes,  to  be  involved  in  one 
observation  equation.  But  this  is  only  a  part  of  the  task,  and  an  easy  part 
at  that;  the  number  of  P.M.  in  such  an  equation  can  be  made  almost  arbitrar¬ 
ily  small  by  decreasing  the  dimension  of  a  spherical  cap  centered  on  an 
observation  point  beyond  which  the  P.M.  are  ignored.  Clearly,  this  repre¬ 
sents  a  simplification  which  should  not  be  overused  because  the  rigor  of 
the  solution  could  suffer.  Besides,  if  the  numbering  of  the  parameters  were 
done  in  a  haphazard  way,  the  span  between  the  P,M.  with  the  lowest  number 
and  the  one  with  the  highest  number  in  the  observation  equation  could  still 
be  great  and  thus  the  bandwidth  would  be  large.  There  might  be  a  great 
number  of  zeros  within  the  band,  but  this  represents  no  advantage  since  such 
zeros  cannot  benefit  the  algorithm.  The  main  task,  therefore,  is  that  of  a 
proper  design  of  the  P.M.  network  and  of  P.M.  numbering. 
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The  most  natural  distribution  of  P.M.  is  that  of  an  equilateral 
grid,  in  that  it  does  not  favor  any  particular  region  and  the  resolution  is 
uniform.  A  basic  example  of  an  equilateral  grid  of  P.M.  appears  below.  The 
parameters  are  numbered  in  columns  (N-S  direction)  because  it  is  anticipated 
that  the  longitudinal  extent  of  the  grid  will  be  greater  than  its  extent  in 
latitude  and  thus  a  smaller  span  will  be  involved  with  one  observation  point. 

If  the  ocean  basin  is  too  great  for  a  strip-like  grid  of  P.M.  to  cover  it, 
overlapping  grids  can  eventually  result  in  a  uniform  geoidal  representation 
of  this  basin  and,  indeed,  of  all  of  the  world's  oceans.  The  observations 
for  each  strip  would  have  to  overlap  as  well,  for  the  sake  of  a  smooth 
transition  of  geoid  contours,  etc.,  between  the  strips.  As  a  practical  matter, 
several  adjacent  strips  may  be  present  in  one  adjustment.  In  a  number  of 
cases  along  the  borderlines,  two  P.M.  parameters  will  be  assigned  the  same 
physical  location,  one  P.M.  belonging  to  the  upper  strip  and  the  other  one, 
with  a  much  higher  rank  number,  belonging  to  the  lower  strip.  Because  of  the 
rank  disparity  the  algorithm  will  treat  these  double  P.M.  completely  inde¬ 
pendently.  To  assure  continuity  between  the  strips  (predictions,  contour 
lines),  each  such  location  could  be  assigned  one  final  parameter  represented 
by  the  average  of  the  two  P.M.  magnitudes.  However,  this  would  imply  the 
renumbering  of  the  parameters  before  the  predictions  could  begin.  The 
simplest  solution  to  this  problem, yielding  exactly  the  same  results,  is  to 
leave  the  numbering  unchanged  and  to  divide  the  magnitude  of  each  double 


1  n^+1  2n^+l  3n^+l 

X  X  X  X  X  X  X 

2  n^+2  2n^+2 

X  X  ©  X  X  X 

3  n^+3  2nj+3 

X  X  X  X  X 

4 

X  X  X  X  X 


2n^  3n^ 

XXX 


In  the  illustration  above  the  "x"  represent  the  P.M.  and  "©"  is  the 
observation  point  whose  observation  equation  is  being  formed.  If  s  is  the 
separation  of  neighboring  P.M.  and  if  the  radius  around  ©  beyond  which  the 
P.M.  are  ignored  in  this  equation  is  just  under  1.5s,  the  span  of  the  P.M. 
involved  is  2n1+3,  the  present  case  comprising  the  P.M.  numbered  1  as  the 
lowest  number  and  Zn^+3  as  the  highest  number.  There  are,  of  course,  many 
P.M.  in  between  which  are  omitted  from  this  observation  equation  (e.g.  those 
numbered  4,  5,  ...,  n1+4,  etc.).  The  number  of  "rows"  of  P.M.  is  denoted  as 
nx  and  the  number  of  "columns"  is  denoted  as  n2.  The  total  number  of  para¬ 
meters  (n)  in  such  a  grid  is  then 

n  =  nxn2  .  (3.1) 


On  the  other  hand,  the  bandwidth  (b)  in  the  matrix  of  normal  equations  under 
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the  circumstances  illustrated  above  is 


b  =  2n.+3  .  (3.2) 

X 

The  bandwidth  is  defined  as  the  largest  number  of  elements  in  any  row  of  N, 
starting  with  the  diagonal  element,  beyond  which  all  the  elements  are  zeros. 
It  is  given  by  the  largest  span  in  the  matrix  of  observation  equations  (A) 
which,  in  the  present  case,  is  2njL+3. 

The  physical  limitations  to  the  geoidal  representation,  both  in 
the  extent  of  the  area  and  in  the  resolution  capabilities,  are  given  by  the 
number  of  elements  in  N  which  have  to  be  stored.  In  the  case  of  no  band¬ 
width  implementation  this  storage  space  (S)  is 

S  =  n(n+l)/2  (3.3) 

elements,  provided  only  the  upper  triangular  portion  of  N  is  stored  (for  the 
whole  matrix  N  this  storage  would  be  n  ).  On  the  other  hand,  if  only  the 
bandwidth  is  stored,  starting  with  the  diagonal  elements,  the  storage  space 
(S')  becomes 

S'  =  bn  .  (3.4) 

This  storage  is  actually  conservative,  in  the  sense  that  the  band  has  been 
extended  by  (b- 1 ) b/2  elements  (zeros)  in  order  to  conform  to  a  rectangular 
form  facilitating  the  implementation  of  a  modified  Choleski  algorithm. 


In  considering  (3. 1) -(3.4) ,  the  storage  savings  associated  with 
the  bandwidth  approach  can  be  expressed  as 

S'/S  =  4(n  +  3/2)/(n+  1)  =  4/n2  .  (3.5) 

The  approximate  relation  in  (3.5)  can  be  adopted  for  a  reasonably  large  nx , 
perhaps  20  or  more.  It  shows  that  with  an  increasing  number  of  "columns"  of 
P.M.  the  storage  economies  grow.  Thus  if  n2 =  100,  the  economies  should  im¬ 
prove  about  25-fold.  This  can  be  illustrated  with  an  example  where  n1=20, 
n^=100  and  thus  n=2000.  The  approximate  value  of  S'/S  was  just  shown  to  be 
1/25  while  the  exact  value  is  1/23.3.  In  particular,  5=  2,001,000  and 
S'  =86,000.  Thus  one  can  see  the  beneficial  effect  of  the  bandwidth  approach 
not  only  from  the  economy  point  of  view  in  itself  but,  especially,  because 
it  allows  an  adjustment  on  a  large  scale  which  would  otherwise  be  unthinkable. 

Due  to  possible  irregularities  in  the  P.M.  grid  the  band  can  ini¬ 
tially  be  larger  than  the  final  b.  If  the  former  is  adopted  in  the  algorithm, 
one  or  more  of  its  first  rows  could  contain  only  zeros.  This  would  detract 
from  possible  economies  but,  most  of  all,  it  would  make  the  algorithm  break 
down,  as  if  N  were  singular.  For  this  reason  it  is  useful  to  search  for  zero 
rows  in  the  band,  eliminate  them  and  arrive  at  a  final  bandwidth.  The 
original  bandwidth  can  be  denoted  as  max  b  and  it  can  be  stipulated  as  2ni+3 
or  perhaps  even  higher.  The  final  bandwidth  is  denoted  as  b;  in  most  cases 
it  will  be  the  same  as  max  b,  but  its  implementation  is  a  simple  matter  well 
worth  the  effort. 

Next,  the  formation  of  the  special  matrix  T  will  be  addressed,  based 
on  one  row  (Ai-row)  of  A.  The  transposition  of  this  row  results  in  a  column 
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which,  when  post-multiplied  by  Ai~row  fills  the  matrix  N  (mostly  with  zeros) 
When  this  operation  is  performed  with  all  the  other  rows  corresponding  to 
all  the  other  observations,  where  the  contributions  are  added  algebraically, 

the  final  matrix  N  is  formed.  All  the  A. -rows  are  assumed  to  have  been 

1 

scaled  by  the  a  priori  sigma  of  the  observations  considered  mutually  inde¬ 
pendent  so  that  weighting  need  not  be  considered.  The  right-hand  side 
(column  vector  U)  of  normal  equations  is  formed  in  a  similar  fashion,  except 
that  the  transposed  A  -row  is  post-multiplied  by  the  corresponding  constant 
term  c  (an  element  of  Lb-L°  in  adjustment  notations).  Due  to  the  span  con¬ 
siderations,  the  Ai~row  will  have  only  b  elements  (actually,  one  would  start 
with  max  b  elements  and,  at  the  end,  "squeeze  out"  the  zero  rows  from  T,  if 
any;  but  it  is  not  difficult  to  imagine  the  following  scheme  apply  with  max 


b  replacing  b) , 


A. -row 

i 


const,  term 


1  2  3 
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Depicted  above  is  the  formation  of  N  with  the  aid  of  one  A^row; 
the  ether  such  rows  would  perform  a  similar  role  as  already  explained.  The 
U-vector  is  also  illustrated.  The  zig-zagging  line  indicates  the  upper- 
criangular  portion  of  N  which  is  of  interest.  The  first  element  in  A^-row 
corresponds  to  the  column  i  in  N.  In  the  computer,  the  rank  number  of  the 
P.M.  i  is  retained  for  addressing  in  N  and  the  Arrow's  addressing  begins 
with  1.  It  ends  with  b  which,  in  N,  corresponds  to  i+b',  where 

b ‘  =  b  -  1  .  (3.6) 

The  important  property  of  the  algorithm  being  developed  is  that  the  T  matrix 
will  not  comprise  the  whole  upper  triangular  portion  of  N,  but  only  its  b-band 
(the  zeros  outside  the  band  are  simply  disregarded).  To  facilitate  the  ad¬ 
dressing,  T  is  enlarged  to  include  the  (nonexisting)  elements  in  N  which 
would  fill  a  triangular  matrix  extending  from  the  upper-left  corner  of  the 
latter.  Such  a  triangular  matrix  would  continue  the  border  of  N  along  the 
dashed  line  for  the  length  of  b'  elements;  its  base  would  also  have  b'  ele¬ 
ments  so  that  the  total  number  of  elements  in  it  would  be  (b-l)b/?..  If 
this  matrix  is  adjoined  to  the  T-portion  of  N,  the  resulting  array  T  will 
have  b  "jagged"  rows  and  n  columns.  This  array  is  finally  "straightened  out" 
to  form  a  rectangular  array  which  is  the  desired  matrix  T  with  b  rows  and  n 
col umns . 

The  rows  in  N  are  represented  by  ascending  diagonals  in  T  as  in¬ 
dicated  below,  where  the  elements  "x"  correspond  to  the  elements  so  denoted 
in  the  last  illustration.  The  ( b-1 ) b/2  nonexisting  elements  of  N  discussed 


above  are  represented  by  zeros  in  T.  Clearly,  the  diagonal  elements  in  N 
occupy  the  last  row  in  T.  Thus  the  "ii"  element  in  N  corresponds  to  the 
"bi"  element  in  T  as  illustrated. 


123...  b'  i  ...  i+b'  ... 


The  matrix  T  is  not  formed  from  N  (this  would  entail  a  large  storage) 
but  directly  from  the  Arrows.  The  element  ii  in  N  has  the  location  bi  in  T; 
the  element  next  to  the  right  in  N,  i,(i+l),  has  the  location  (b-l),(i+l)  in 
T,  followed  by  (b-2),(i+2),  etc.  The  last  element  within  the  band  b  in  that 
row  of  N  has  the  location  l,(i+b')  in  T.  Thus  the  matrix  N  has  served  only 
mentally,  without  any  computer  core  requirements.  All  the  "ascending  diagonals 
in  T  are  formed  according  to  the  same  principle,  except  that  the  last  diagonals 
end  at  the  column  n  and  do  not  reach  the  first  row  of  T  (they  have  fewer 
elements  than  b).  The  vector  U  is  not  modified  in  conjunction  with  the  for¬ 
mation  of  T.  Mathematically,  the  above  principle  of  addressing  can  be 
symbolized  as  follows: 
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(3.7a) 


i  ,i  b  *  i  *  i 

i  ,j  b-(  j-i )  ,j ;  j  =  i+1  ,i+2 . n.  (3.7b) 

Next,  the  modified  Choleski  algorithm  will  be  developed  to  treat 
the  problem  at  hand,  in  conjunction  with  the  above  matrix  T.  With  the  custo¬ 
mary  adjustment  notations  the  least-squares  solution  of  the  vector  of  n 
(unknown)  parameters,  X,  is  expressed  as 

X  =  N— 1U  , 

where 

N  =  ATA  ,  U  -  AT(Lb-L°)  , 

and  where  Lb  is  the  vector  of  observations  while  L°  is  the  vector  of  the  cor¬ 
responding  initial  values.  As  mentioned  earlier,  the  weighting  of  observa¬ 
tions  need  not  be  applied  here.  The  variance-covariance  matrix  (E)  for  the 
parameters  after  the  adjustment  is  obtained  as 

Ex  -  N-1  . 

If  F  is  a  w-vector  of  linear  combinations  of  the  parameters, 

F  =  iFx  , 

where  U  is  a  matrix  of  dimensions  n>w  ,  one  has 

t  =  IFn-1!!  , 

V 
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which  follows  from  the  law  of  the  variance-covariance  propagation. 

If  N  is  nonsingular  it  is  positive  definite  and  it  can  be  expres¬ 
sed  as  a  product  of  a  nonsingular  matrix  by  its  transpose,  namely 

N  =  TTT  • 

In  the  "regular"  Choleski  algorithm  (without  the  bandwidth  considerations) 

T  is  in  general  an  upper- triangular  matrix  with  no  other  particular  pattern 
and  thus  no  special  system  of  addressing.  But  this  does  not  alter  the 
present  general  considerations.  In  taking  advantage  of  the  triangular  de¬ 
composition  the  following  results  are  obtained  at  once: 

X  =  T-1R  ,  r  =  R1^  , 

r 

where 

R  =  (Tt)“1U  ,  R  »  (t1)-1!)  , 

and  where  the  property  N-1 =  T-1(Tt)-1  has  been  utilized. 

In  collecting  the  pertinent  formulas,  one  first  writes 

TtT=  n  ,  ttr=  u  ,  TTR  =  u  , 

which  can  be  represented  by 

TX[  T  I  R  •  R  ]  =  [  N  j  U  j  U  ]  ,  (3.8a) 

n«n  n«l  n«w  n»n  n-1  n.w 

where  the  dimensions  are  also  indicated.  The  remaining  formulas  of  interest 
are 
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X  =  T"1  R  ,  (3.8b) 

n»l  n«n  n*l 


w.w  w.n  n«w 


The  computation  of  the  elements  in  T,  R  and  R  via  (3.8a)  is  sometimes  called 
"forward  reduction"  and  the  computation  of  X  in  (3.8b)  may  be  called  "back 
substitution". 

The  elements  of  T  in  (3.8a)  are  computed  row  by  row,  starting  with 
the  "ii"  element  and  proceeding  to  the  right  until  the  "in"  element  is  com¬ 
puted.  The  computation  of  the  elements  in  R  and  R  follows  the  same  alge¬ 
braic  pattern,  to  the  extent  that  R  and  the  columns  in  R  can  be  treated 
simply  as  additional  columns  to  T  (beyond  the  n-th  column).  If  the  computa¬ 
tion  proceeds  as  indicated,  i.e.,  starting  with  the  first  element  in  the 
first  row  and  following  through  with  that  row  till  completion,  then  start¬ 
ing  with  the  "22"  element,  etc.,  the  elements  in  T,  R  and  R  thus  computed 
can  replace  the  corresponding  elements  in  N,  U  and  U  because  the  original 
elements  in  the  latter  three  arrays  will  no  longer  be  needed  in  the  computa¬ 
tions.  This  leads  to  savings  of  the  core  space  where  only  N  (upper  tri¬ 
angular  portion),  u  and  U  need  to  be  stored,  becoming  later  T  (upper- 
triangular),  R  and  FL  During  the  step  (3.8b),  the  elements  of  X  may 
similarly  replace  the  elements  of  R  and  thus  the  column  which  started 
originally  as  U  will  become  X  at  the  end  of  all  the  computations.  The 
Choleski  algorithm  is  then  seen  as  performing  the  following  transformations 
of  elements  in  the  same  computer  storage: 


-38- 


[N  U  U]-  [T  R  R]  -*■  [ T  X  R]  . 


(3.9) 


The  algebraic  operations  carrying  out  (3.8a)  in  the  "regular" 
Choleski  algorithm  (with  "regular"  addressing  in  T)  are 


t. . 
11 

=  (n..  -  t2  -  t2  -  ...  -  t2  .)*  , 

11  li  2i  i-l,i 

(3.10a) 

t. . 

=  (n. .  -  t. .t. .  -  t_.t_.  -  ...  -  t.  ,t.  ,  . )/t. .  , 

lj  li  lj  2i  2j  i-l,i  1-1,3  n 

(3.10b) 

r. 

i 

=  (u. -  t.  .r,  -  t_.r_-  ...  -  t.  ,  .r.  ,)/t. .  , 

1  li  1  2i  2  1-1,1  1-1'  n 

(3.10c) 

Fij' 

,=  (IT.  .  -  t.  .7.  . ,  -  t_.r_  t.  .  .7.  .  . ,  )/t. .  . 

1  13 •  li  13'  2i  23'  i-l,i  1-1,3'  li 

(3.10d) 

Here  j  starts  with  i+1  and  proceeds  to  n  whereas  j'  proceeds  from  1  to  w, 
all  for  the  given  i;  subsequently  i  increases  by  one,  and  so  on  (i  started 
with  1  and  will  end  with  n).  In  the  step  (3.8b),  i  starts  with  n,  then 
continues  with  n-1,  etc.,  and  ends  with  1.  Thus  x  is  computed  first,  x  , 

n  n-i 

follows  in  the  process  where  the  known  x  is  also  used,  etc.,  hence  the 

n 

name  "back  substitution".  With  this  provision  (3.8b)  translates  into 


x. 

1 


(r.  -  t.  x  -  t.  .x  ■ 
i  m  n  i,n-l  n-1 


i+lXi+l 


which  can  start  after  the  completion  of  the  step  (3.8a). 


(3.11) 


The  next  step  is  only  intermediate  and  corresponds  to  (3.10a-d) 
and  (3.11)  in  the  case  of  a  banded  N  matrix,  but  without  any  change  in 
addressing.  Advantage  is  merely  taken  of  the  zeros  outside  the  band  in 
that  they  are  disregarded  and  not  used  in  computations.  A  remark  to  be 
made  with  regard  to  the  indices  is  that  in  (3.12a-d)  below  they  can  be 
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only  positive  (zero  and  negative  indices  are  disregarded),  that  j  in 
(3.12b)  ends  with  i+b'  or  n,  whichever  is  smaller,  and  that  indices  larger 
than  n  in  (3.13)  are  disregarded. 


2  2 

t. .  =  (n. .  -  t.  .  ,  .  -  t.  .  _  .  -  . 

ii  ii  i-b+1,1  i-b+2,i 


-  t2  )  2 

i-l,i/  * 


(3.12a) 


t.  .  =  (n.  .  -  t_  .  .t.  .  .  .  -  t.  .  _  .t.  __  t.  .  ,t.  .  ,)/t. .  ,  (3.12b) 

l]  13  3~b+l,i  3~b+l,3  3~b+2,i  3~b+2,3  i-l,i  1-1,3  11 


r.  =  (u.-t.  ,  .  r.  _  ,-t.  ,  _  .  r.  .  -  ...  -  t.  _  .r.  _)/t. .  , 

i  i  i-b+1,1  i-b+1  i-b+2,i  i-b+2  i-l,i  l-l  n 


(3.12c) 


rij‘  ^  Ui  j  *  ^i-b+l,iTi-b+l, j '  ^i-b+2 ,iri-b+2,  j  '  ***  j  ’  ^^ii’ 


X.  =  (r.  -  t.  .  ,X.  ,  -  t.  . , ~X. . -  -  . . .  -  t.  .X.^.  ,)/t. .  . 

i  i  i,i+l  i+l  i,i+2  i+2  i,i+b-l  i+b-1  ii 


(3.12d) 

(3.13) 


Finally,  advantage  is  taken  of  the  banded  structure  of  N  by  replacing 
it  with  T  described  previously.  This  entails  the  change  in  addressing  pre¬ 
sented  as  (3.7a,b).  The  matrix  T  then  also  changes  its  shape  and  addressing; 
for  this  reason  it  should  perhaps  be  called  by  a  different  name,  but  here  the 
symbol  "T"  is  retained  with  the  knowledge  that  it  is  now  a  rectangular  matrix 
of  dimensions  b«n.  This  having  been  said,  the  elements  in  the  (new)  T  matrix, 
in  R  and  R,  and  finally  in  X  are  presented  below  where,  in  addition  to  (3.6), 
also  the  following  notation  is  used: 

b"  =  b'  -  (j  -  i)  .  (3.14) 

The  formulas  are  written  in  their  final  form  convenient  for  programming.  The 
principal  index  i  runs  again  from  1  to  n,  the  index  j  within  each  i  ends 
again  with  inf(i-b',n),  etc. 


-40- 


3.15a) 


b- (3-1)  ,3  9“.  :-!+£, i  ay  bi 


,3.15b) 


i  <b :  r . 


1  ”  ^Ui  ’  tb-i+4,ir£^tbi  ’ 


(3.15c) 


l>b:  r4  =  (u.  -  ^  tur..b+t)/tbl  . 


(3.15c') 


i<b:  r...  =  (u...  -  T  t.  .  .  .r„.,)/t  .  , 

13*  13 '  b-i+£,i  S.2  bi 


( 3 . 15d) 


i  >b :  r . 


ij'  ^ Ui j  •  "  ^  tHri-b+£,j,^tbi 


( 3 . 15d ' ) 


inf (b* ,n-i) 


xi  ■  -  I 


^b-£ ,  i+£Xi+£  ^  ^bi  * 


(3.16) 


Instead  of  the  full  variance-covariance  matrix  of  F,  only  the  chosen 
sigmas  and  the  correlation  coefficients  will  be  computed.  The  input  para¬ 
meters  will  be  the  pairs  jj  and  which  will  serve  to  choose  the  appropriate 
columns  of  If  and  compute  the  desired  values  as  follows: 


/  r  -2  A 

=  (  l  r  )2  , 

31  £=1 


or-(l  r2  )H  , 

J 2  £=1  J2 


Pj,i,=  (  l  r..,r9.,)/(a.,aA,) 


(3.17a) 


(3.17b) 


(3.17c) 


Some  economies  can  be  realized  if  one  takes  into  account  that  several  elements 
in  the  columns  of  U  may  be  zeros;  however,  only  those  zero  elements  which 
precede  the  first  nonzero  element  are  beneficial.  Similar  considerations 
could  apply  also  with  regard  to  U  but  this  is  merely  a  theoretical  case 
since  as  a  rule  U  is  a  full  vector.  Suppose  that  in  a  column  j1  of  TJ  the 
first  nonzero  element  is  at  the  address  i=J‘.  Then  the  first  nonzero  element 
in  the  j'-th  column  of  R  appears  also  at  the  address  i=J'  so  that  the  com¬ 
putations  of  the  preceding  elements  as  in  (3.15d),  ( 3 . 15d ' )  can  be  skipped 
and  these  elements  can  be  set  directly  to  zero.  Furthermore.the  multiplica¬ 
tions  on  the  right-hand  sides  of  (3.15d),  ( 3 . 15d ' )  with  these  zero  elements 
can  likewise  be  skipped,  in  other  words,  2,=1  in  the  summations  can  be  re¬ 
placed  with  the  appropriate  J'.  In  the  same  vein,  «=1  in  (3.17a,b)  can  be 
replaced  by  the  corresponding  and  0^  »  respectively,  and  «,=1  in  (3.17c)  can 

be  replaced  by  supfJ^J^). 

Finally,  one  notices  that  the  symbolism  (3.9)  would,  in  the  modified 
Choleski  algorithm,  read  as 

TUU  +  TRR  +  TXR.  (3.18) 

b«n  n«l  n«w  b«n  n«l  n«w  b.n  n«l  n«w 

The  brackets  have  been  left  out  because  the  complete  expression  is  not  one 
matrix. 

The  approach  utilizing  the  P.M.  parameters  is  based  on  the  residuals 
from  a  previous  global  adjustment  in  terms  of  spherical -harmonic  potential 
coefficients.  The  P.M.  approach  considers  the  residuals  in  the  exact  loca¬ 
tions  on  the  globe  where  the  actual  altimeter  measurements  took  place. 
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There  are  other  interpolation  and  approximation  techniques  where  the  obser¬ 
vations  are  "translocated"  to  some  strategically  advantageous  positions, 
such  as  to  grid  locations  (equilateral,  geographical,  etc.).  Such  techniques 
may  be  very  crude  (the  location  is  considered  shifted  but  the  value  is  taken 
without  any  correction)  or  quite  elaborate  (e.g.,  as  in  the  "best  linear 
prediction"  method).  Since  the  values  at  shifted  locations  are  estimated 
and  not  measured,  any  such  technique  entails  approximations  which  are  avoided 
in  the  P.M.  adjustment.  As  a  trade-off,  some  of  the  above  techniques  are 
very  efficient  computationally.  The  present  approach  with  the  modified 
Choleski  algorithm  seeks  to  improve  the  computational  efficiency  while  re¬ 
laxing  the  rigor  of  the  solution  as  little  as  possible. 


4.  MODEL  FOR  DEFLECTION  RATES  IN  TERMS  OF 
SPHERICAL  HARMONICS  AND  POINT  MASSES 

The  need  to  develop  mathematical  foundations  for  new  observational 
modes  is  dictated  by  the  surge  in  various  new  or  improved  operational  systems. 
For  example,  broad  advances  have  been  made  in  the  area  of  inertial  instrumen¬ 
tation  which  has  been  used  with  growing  success  in  inertial  navigation.  In 
using  this  instrumentation,  the  observational  process  can  be  designed  to  yield 
differential  changes  in  the  deflections  of  the  vertical  along  a  given  route. 
Thus  the  quantities  to  be  modeled  by  a  given  set  of  parameters,  such  as  the 
spherical -harmonic  (S.H.)  potential  coefficients  or  the  point  mass  (P.M.) 
magnitudes,  can  be  viewed  as  the  deflection  rates  at  chosen  (observational) 
points.  Such  points  would  usually  be  distributed  along  profiles  whose  azi¬ 
muth  is  denoted  by  a.  At  each  point  the  measured  quantities  to  be  modeled  are 
denoted  as  £  and  rj,  associated  with  a  specific  a.  This  type  of  data  involves 
second  derivatives  of  the  disturbing  potential  (T)  with  respect  to  length 
elements  (ds),  which  could  be  particularly  useful  for  the  resolution  of 
short-wavelength  features  in  a  local  gravity  field.  The  development  contain¬ 
ed  herein  is  based  on  [Hotine,  1969]  abbreviated  as  [H]  and  on  [Blaha,  1980] 
abbreviated  as  [B] .  It  could  be  considered,  in  fact,  to  be  an  extension  of 
Chapter  4  in  [b]  insofar  as  the  S.H.  and  P.M.  (single  layer)  parameters  are 
concerned,  and  it  could  follow  Section  4.6  in  this  reference. 


4.1  Mathematical  Background 


Similar  to  the  material  contained  in  paragraph  19  on  page  7  of 
[H] ,  the  derivative  of  a  function  F  with  respect  to  the  length  element  along 
a  desired  direction  reads 

3F/3S  =  (3F/3xr)(3xr/3s)  =  Frpr  , 

which  is  an  invariant  (independent  of  the  coordinate  system  represented  by 
xr) ,  where 


=  gradient  of  F  , 

pr  =  unit  vector  in  the  desired  direction, 
whose  length  element  is  ds. 

In  the  present  context,  pr  lies  in  a  level  surface  and  its  azimuth  is  a. 
Similar  to  equation  12.006  of  [h]  we  write 

pr  =  A'  si  not  +  Zr(2)  cosa  , 

where  (1)  refers  to  an  easterly  direction  along  which  the  differential  dis¬ 
tance  is  denoted  ds  and  (2)  refers  to  a  northerly  direction  associated  with 
ds.,.  In  the  spherical  coordinate  system  {A,*,r},  the  first  two  coordinates 
are  the  geocentric  longitude  and  latitude,  respectively,  and  the  third  is 
the  geocentric  radial  distance.  In  agreement  with  (3.48)  of  [8]  we  have 

Aj1}  =  { 1/ ( r  cos*),  0,  0},  Zr{2)  =  (0,  1/r,  0}, 

pr  =  {(1/r  cos*)sina,  (l/r)cosa,  0}  ,  (4.1) 
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where  the  index  "r"  should  not  be  confused  with  the  radial  distance  denoted 
by  the  same  letter.  Here  and  in  the  sequel,  (1/r  cos  o') ,  etc.,  are  written 
to  represent  [l/(r  coso)] ,  etc.,  since  the  latter  notation  would  sometimes 
contribute  to  a  profusion  of  parentheses  and  brackets;  there  is  no  need  for 
the  former  notation  to  be  confused  with  ( l/r)cos<jj".  As  is  customary  with 
tensor  notations,  the  summation  convention  applies. 

If  3F/3s  represents  a  deflection  of  the  vertical  (the  sign  change 
follows  from  the  convention  explained  on  page  79  of  [B]),  F  becomes  the  geoid 
undulation  N  and  we  write  for  the  deflections  in  the  north  and  east  direc- 


tions,  respectively: 

S  =  -3N/3s  =  -N  A'  i  =  '(l/r)SN/3$  , 

N  r  (Z) 

(4.2) 

n  2  -SN/Ss^.  =  -N  l*.  =  -{1/r  cos^)3N/3X  . 

E  r*  v 

(4.3) 

The  rates  of  these  quantities  along  the  general  azimuth  a,  along  a=0  and 
along  a-  90°  are  expressed  as 


k  =  3£/3s  =  -32N/3s3s  =  i  pr  , 


N 


n  =  3n/3s  =  -32n/3s_3s  =  n  Pr  ; 

ht  JT 


(4.4) 

(4.5) 


L  =  3F/3sm  2  -32N/3S2  =  ?  ir 

N  N  N  ( 


(2) 


\  ;  ;  -S'N/9se3sh  -  nr^2); 


E 


v;/3sr  =  -32N/3sn3se  -  , 


rE  *  3^3SE  "  ^  N/3SE  =  V(l) 
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The  above  relations  do  not  mean  that  the  same  formulas  will  be 


recovered  for  t  and  n.  Instead,  another  relationship  will  be  derived  for 
checking  purposes.  First  suppose  that  scalar  functions  of  position  F  and  H 
are  such  that 

H  =  3F/3s  =  F  pr  . 

r 

Next  identify  another  direction  and  its  unit  vector,  etc.,  by  primes,  namely 
H'  =  3F/3s 1  =  F  p,r  . 

r 

Two  kinds  of  derivatives  are  now  formed: 

3H/3S '  =  32F/3s3s 1  =  (F  pr)  p,s  , 

r  s 

3H73S  =  32F/3s  ' 3s  =  (F  p,r)  ps  . 

r  s 

Using  the  rules  for  covariant  differentiation  according  to  Chapter  3  of  [H] , 
we  have 


(F  pr)  =  F  pr  +  F  pr  , 
r r  s  rsr  r  s  * 


Ps  ■  3‘>r/s’'s  +  r>  • 


the  symbols  Frg,  Tgk,  symmetric  in  the  lower  indices,  eventually  cancel 
out  and  need  not  be  elaborated  ?on.  Straightforward  algebra  yields 


32F/3s'3s  -  32F/3s3s 1  =  Fr[(3p,r/3xs)ps  -  ( 3pr/3xs)p ,s] . 


(4.6) 


We  now  associate  F  with  N  and,  anticipating  the  spherical 


approximation,  write 


N  =  ( 1/G ) T  , 


(4.7) 


where  6  is  the  average  value  of  gravity  at  the  earth's  surface  (approximately 
980  gals).  Equation  (4.7)  follows  from  the  Bruns  formula  presented  on  page  85 
of  [Heiskanen  and  Moritz,  1967].  If  ds  and  ds ‘  are  identified  with  dsN  and 
ds  .  respectively,  and  other  quantities  are  identified  accordingly,  it  is 
readily  deduced  that 

32N/3s  3s  -  32N/ 3s  3s  =  L  -  n  =  ( l/Gr2COS?)  tg^  3T/3A  .  (4.8) 

E  N  N  E  E  N 


In  analogy  to  G  in  (4.7),  r  in  the  final  formula  can  be  substituted  for  by  R, 
the  earth's  mean  radius  (approximately  6371  km),  for  the  deflection  rates  as 
modeled  customarily  for  the  earth's  surface. 

In  considering  (4.1)  and  (4.7),  the  formulas  (4.2)  -  (4.5)  yield 


5  =  -(1/Gr2)  [(l/cos$>)(32T/3A34>)sina  +  (32T/ 34>2)cosa]  , 


2t  /  iX2  ’ 


(4.9) 


f,  =  -(l/Gr2cos4>)[(l/coscJ>)(32T/3A2)sina  +  (tg<p  3T/3X 


A  \  /  *>  \  2  > 


+  32T/3A3<j>)cosa].  (4.10) 

However,  if  the  computation  of  32T/3^2  is  sought  to  be  avoided  as  in  the  S.H. 
model,  one  can  take  advantage  of  the  Laplace  condition  for  harmonic  functions, 
in  particular, 


AT  =  0  . 
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The  alternate  formula  for  i  then  reads 

(4. 9') 

k  =  - ( 1/G) { ( 1/ r  cos0)(32T/3A3$)sinct  +  [-(l/r2cos2^)32T/3X2  -  32T/3r2 

+  (tg?/r2)3T/3$  -  (2/r)3T/3r]cosa)  , 

where  the  expression  inside  the  brackets  replaces  32T/3^2  in  (4.9).  Upon 
setting  a=90°  in  (4.9)  or  (4.9'),  and  a=0  in  (4.10),  the  verification  equation 
(4.8)  follows.  If  the  formulas  for  the  deflection  rates  in  spherical  approxi¬ 
mation  should  be  used  at  aircraft  or  satellite  altitudes,  r  and  G  values  above 
would  be  substituted  for  by  their  counterparts  for  the  geop  in  question. 


4.2  Spherical -Harmonic  Approach 


In  spherical  approximation,  the  selected  formulas  in  (3.77)  of 
[B]  for  r=R  can  be  rewritten  in  terms  of  a  set  of  S.H.  potential  coefficients 
truncated  at  n=N  as 


3T/3A  =  RG  l  AS'(n)  , 

n=2 


N  _ 

3T/3<JT  =  RG  l  AS(n)  , 

n=2 


3T/3r  =  -G  l  (n+l)AS(n)  , 

n=2 


32T/3A2  =  RG  l  AS"(n)  , 

n=2 


32T/3A3^  =  RG  l  AS'(n)  , 

n=2 


32T/3r2  =  (G/R)  £  (n+l)(n+2)AS(n)  , 

n=2 


(4.11) 


where  the  notations,  according  to  (3.76)  and  (3.78a-d)  of  [B] ,  are 


AS(n)  =  7  (AC  cos  mA  +  AS  sin  mA)P  (sinljf)  , 

n  nm  nm  nm 

m=0 


AS'(n)  =  3AS/3A  =  l  m(-AC  sin  mA  +  AS  cos  mA)P  ( s i n <{> )  , 

u  r,  nm  nm  '  nm  T 

m=u 
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(4.12 


AS"(n)  =  9  AS/9A  =  -  £  m2(AC  cos  m\  +  AS  sin  mX)P  (sin<}>)  , 

_  „  nm  nm  nm 

m=u 


AS(n)  =  9AS/9<F  =  l  (AC  cos  mX  +  AS  sin  mX)dP  (sin^)/d<jr  , 

nm  nm  nm  Y  * 

m=u 


-  2  _  n  _  _ 

AS'(n)  =  3  AS/9X3cj>  =  Y  m(-AC  sin  mX  +  AS  cos  mX)dP  (sin<)>)/d<t> 

_  nm  nm  nm 

m=0 


With  C*q  representing  the  reference  field,  the  notations  are 


AC  =  C  -  C* 

no  no  no 


for  n=2,  4  and  perhaps  6,  the  remaining  AC  and  all  of  the  AS  being  the 

nm  nm 

coefficients  C  and  S  themselves. 

nm  nm 

Upon  utilizing  the  notations  (4.11),  (4.12),  the  formulas  (4.9') 
and  (4.10)  lead  to 

k  -  -( 1/R) { ( l/cos<p)  l  AS'(n)  sina  -  [(1/cos  <F)  £  AS"(n) 


N  _  N  _ 

+  l  n(n+l)AS(n)  -  tg<j>  l  AS(n)]  cosa)  , 

n=2  n=2 


N  N 

n  =  -(1/R  cos^){(l/cos^)  l  AS"(n)  sina  +  [tg^  l  AS'(n) 

n=2  n=2 

N  _ 

+  l  AS'(n)]  cosa)  . 

n=2 


(4.13) 


(4.14) 


If  a=90°  and  a=0  are  used  in  (4.13)  and  (4.14),  respectively,  the  verification 
equation  (4.8)  follows  upon  utilizing  the  notation  from  the  first  equation 
of  (4.11). 
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In  following  the  same  sequence  of  derivations  as  in  [B] ,  the 
partial  derivatives  of  the  above  quantities  with  respect  to  the  S.H.  coef 
ficients  are  expressed  as 

3C/9Cn0  =  (1/R)  [n(n+l)Pn(sin<JT)  -  tg$  dPn(sin^)/d$]cosa  , 

3£/3C  =  ( 1/R) ( A  sin  mX  -  B  cos  mX)  , 

nm 

3£/3S  =  -( 1/R) ( A  cos  mX  +  B  sin  mX)  ; 

nm 

9n/3C  *  0  , 

no 

3n/3C  =  (m/R  cos^)(C  cos  mX  +  D  sin  mX) , 

nm 

3fi/3S  =  (m/R  cos^)(C  sin  mX  -  D  cos  mX)  ; 

nm 

the  new  notations  in  these  expressions  are 
A  =  (m/cos^)[dPnm(sin<)>)/d$]sina  , 

B  =  { [(l/cos2$)m2  -  n(n+l)]P  (sini)  +  tg<f  dP  (sin^)/d^}cosa  ; 

L  J  nm  nm 

C  -  (m/cos<j)’)P  (si n^>")  sina  , 

nm 


D  =  f tgcp  +  dPnm(sin$)/d4>Jcosa  . 


The  notations  in  this  analysis  will  be  the  same  as  those  used  in 
[B] ,  pages  67  and  68  for  the  point  masses  located  at  the  depth  d  below  the 
sphere  of  radius  R: 


In  the  case  the  observation  points  should  not  be  at  the  earth's  surface 
but  at  aircraft  or  satellite  altitudes  (then  R  is  larger  than  6371  km), 
d  may  have  to  be  increased  beyond  the  value  indicated  above  if  the  loca¬ 
tions  of  the  point  masses  are  still  desired  to  be  underneath  the  actual 
earth's  surface. 

The  deflection  rates  at  the  observation  point  i  can  be  transcribed, 
in  spherical  approximation,  from  (4.9)  and  (4.10)  as  follows: 

q  =  -( l/GR2)  [( l/cosd>i)(32Ti/3Xi34>i)sinct  +  (32T./3<*>2)cosa]  ,  (4.16) 

f,  =  -(l/GR2cosq)[(l/cos$.)(32T. /3X2)sina  +  (tg*.  3T./3Xi 

+  32Ti/3Xi34>i)cosa]  .  (4.17) 

The  required  partial  derivatives  of  T  can  be  found  from  (4.15),  and  they 
are  listed  as  (4.7a-i)  in  [8].  Using  these  derivatives,  after  a  few 
algebraic  manipulations  we  find 

q  =  (Rj/GR)  l  (l/^.)[(l/cos<j>i)(pijtij  -  tg<(>i)qijsina 

"  ^ Pi jli j  ’  »  (4.18) 

n.  =  -(R^GR  cosq)  I  (l/ilJj)[(l/cosq)(pijq2j 

-  cos4>. costj). cosAX.  .) si na  -  p.  ,q.  .t.  .  cosal(kM).  ,  (4.19) 

i  j  13  ri]  J  3 


t.  .  *  cos*. sin4> .  -  sin;>.cos<{> .cosAA.  .  , 

13  Y1  YJ  *1  YJ  LJ 

q.  .  =  cos<t>.cos4>.sinAA.  .  . 

3-D  1  j  ij 

Next  we  consider  (4.18)  and  (4.19)  for  one  point  mass  j  and  vary 
the  position  of  the  observation  point  i  by  distances  ns,  where  n  is  the  num¬ 
ber  of  intervals  s  along  a  general  direction.  As  in  [B] ,  this  strategy  makes 
it  possible  to  see  at  how  many  intervals  away  from  the  observation  point  the 
P.M.  parameters  still  significantly  affect  the  modeled  value  (and  vice  versa), 
and  should  be  included  in  the  pertinent  observation  equation  or  prediction 
formula.  The  inclusion  of  all  the  P.M.  parameters  in  a  large-area  adjustment 
would  entail  prohibitive  computer  core  and  run-time  requirements.  The  con¬ 
sideration  of  only  one  point  mass  removes  the  summation  signs  in  (4.18)  and 
(4.19).  As  in  [B]  in  comparable  situations,  the  latitude  of  j  is  chosen  zero, 

<b.=0,  which  here  simplifies  also  the  formulas  for  t.  .,  q.  .  and  cosi(>.  .. 
j  ij  1 j 

Since  we  are  dealing  with  small  changes,  first-order  approximations  are  further 
introduced,  such  as 

cosAA. .  =  1  ,  sinAA.  .  =  AA. .  , 

3-j  30  i} 

cos4>i  =  1  ,  sin4>±  =  A^  , 


etc. ,  where 


In  carrying  out  the  straightforward  algebra  (in  this  process  one 


also  realizes  that  is  about  two  orders  of  magnitude  greater  than  a<j>  ) , 
it  follows  from  (4.18)  and  (4.19)  that 

£.  =  -(R1/6R)(l/£jj)[pijA0>.jAXijsina  +  (p.-A*2.  -  l)cosaJ(kM).  ,  (4.20) 

\  =  -(R1/GR)(l/J^j)[(p..AX2j  -  l)sina  +  p.  .A*.  .AX.  .cosa] (kM) .  .  (4. 20') 


Clearly,  whatever  values  are  obtained  for  5.  with  a  given  azimuth  a=a*  , 

are  also  the  values  obtained  for  f|.  with  a=90°-a'  ,  and  with  A<J>  and 

1  ij 

AX  interchanged.  Therefore,  it  is  sufficient  to  examine  only  one  of 
these  equations,  for  example  (4.20),  with  regard  to  cut-off  considerations 
in  the  "worst  situation".  Such  a  situation  depicts  the  greatest  influence 
of  the  P.M.  at  j  upon  the  modeled  value  at  i  at  various  distances  j -i ,  both 
with  respect  to  a  and  8,  where  0  is  the  azimuth  measured  from  j  to  i. 

In  angular  measure  s  is  denoted  as  Aw  and  the  angular  distances 
between  j  and  i  proceed  in  n-multiples  of  Aw.  We  then  have 


Vi;.  =  n  Aw  ,  A4>ij  =  n  Aw  cos0  ,  AX_  =  n  Aw  sin0  .  (4.21) 


In  view  of  (4.20),  the  following  approximate  relations  are  developed: 


i2.  -  R2Aw2(n2  +  0.64)  , 


pi.  =  3/[Aw2(n2  +  0.64)]  , 


>  (4.22) 


-Sfi- 


where  d=0.8s  from  (4.15')  in  angular  measure  has  been  taken  into  account. 
Equation  (4.20)  thus  becomes 

(4.23) 

=  -C[l/(n2+ 0.64)5/2][3n2sing  cos8  sina  +  (3n2cos2S -  n2 -  0.64)cosa]  , 

C  =  Rx(kM)  ./(GR4Au>3)  .  (4.23') 


In  order  to  find  the  worst  situation,  the  function  f(a,8)  forming 
the  second  brackets  in  (4.23)  is  examined  for  local  extremes,  where 

f(a,6)  =  a  sine  cosB  sina  +  (b  cos2S  -  c)cosa  ,  (4.24) 

a  =  3n2  ,  b  =  3n2  ,  c  =  n2  +  0.64  .  (4.24') 


One  first  forms  3f/3a  =  0  and  3f/3$  =  0  and  seeks  the  solution  for  8,  for  example, 
after  having  eliminated  a.  The  result  is  represented  by 

cos28  =  (2bc  -  a2)/ [2(b2  -  a2)]  . 


Without  any  further  analysis  one  notices  that  such  a  local  extreme  does  not 
exist  due  to  a=  b  in  (4.24').  One  can  then  proceed  by  examining  different 
azimuths  a,  held  fixed,  and  searching  for  an  extreme  with  respect  to  8.  In 
using  this  approach,  the  azimuths  a =0,  45°  and  90°  are  adopted. 

Azimuth  g=0.  Here  the  second  brackets  in  (4.23)  contain  only  the 
second  term  (with  cosa  replaced  by  1).  This  is  then  called  f(8)  whose  local 
extremes  are  sought.  Simple  algebraic  manipulations  reveal  that  an  extreme 
occurs  at  8=0  and  that  it  represents  a  maximum;  the  second  derivative  of 
f(8)  with  respect  to  8  evaluated  with  8=0  is  negative.  However,  we  are 
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interested  in  a  maximum  in  absolute  value,  which  is  seen  to  coincide  with 
the  maximum  just  found  when  f(g)  reaches  at  least  n  +  0.64.  With  3=0 
this  happens  whenever  n>  1.13.  Since  we  are  especially  concerned  with  the 
distances  j-i  in  the  vicinity  of  1.5s,  the  results  implied  by  the  above 
procedure  are  satisfactory.  Thus,  upon  disregarding  the  signs  the  formula 
and  the  values  below  represent  the  absolute  maximum  for  the  intervals  of 
interest: 


a=0,  3=0; 

h  - 

-C[(2n2  -  0.64)/(n2 

+  0.64)5/2 

n  =  0  . 

100% 

n  =  1  . 

,  -20.3% 

n  =  0.5  . 

9.6% 

n  =  1.5  . 

.  -13.9% 

n  =  2  . 

,  -  8.2% 

(4.25) 


l  (4. 25') 


Note:  The  formula  (4.25),  examined  for  n,  reveals  a  primary  maximum  at  n=0 
associated  with  "100%",  and  a  secondary  maximum  at  n=l  (more  precisely 
at  n=0.98) . 


Azimuth  q=45°.  In  this  case  both  sina  and  cosa  in  (4.23)  are 
replaced  by  0.7071  taken  out  of  the  brackets.  A  maximum  for  the  new  f(8) 
occurs  at  3=22.5%  and  it  corresponds  to  the  absolute  maximum  already  for 
n  >  0.89,  which  is  entirely  satisfactory.  Similar  to  (4,25),  (4.251), 
this  configuration  is  represented  by 

a=450  ,  3=22.5°;  =  -0(0.7071)  [(2.62n2  -  0.64ytn2  +  0.64)5/2] ;  (4.26) 
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(4.26') 


n  =  0  100%  n  =  1  -29.5% 

n  =  0.5  -1.0%  n  =  1.5  .  -19.0% 

n  =  0.91  .  -30.1%  n  =  2  -10.9% 

Note:  In  addition  to  the  primary  maximum  at  n=0,  the  formula  (4.26)  has  a 
secondary  maximum  at  n=0.91. 

Azimuth  g=90°.  Now  only  the  first  term  inside  the  second  brackets 
in  (4.23)  remains,  for  which  the  (absolute)  maximum  occurs  at  B=45° .  This 
situation  is  depicted  as  follows: 

a=90° ,  6=45°;  =  -C[l.5n2/(n2+  0.64)5/2]  ; 

n  =  0  .  0  n  =  1  .  -79.8% 

n  =  0.5  .  -  92%  n  =  1.5  .  -43.7% 

n  =  0.65  .  -100%  n  =  2  .  -23.7% 

n  =  2.5 .  -13.8% 

Note:  The  formula  (4.27)  reveals  only  the  primary  maximum  at  n=0.65. 

In  closing  this  part  of  the  analysis,  we  comment  that  due  to  the 
symmetry  in  the  deflection-rate  formulas  in  conjunction  with  the  point  mass 
j  located  in  the  vicinity  of  the  observation  point  i,  it  is  sufficient  to 
consider  the  azimuths  a  and  8  in  the  first  quadrant.  In  examining  three 
representative  azimuths  a  (0,  45°,  90°)  in  conjunction  with  ,  it  is 
noticed  that  a=90°  depicts,  by  far,  the  "worst"  situation.  At  a=0  the  in¬ 
fluence  (at  i)  of  the  point  mass  (at  j)  is  felt  the  least;  the  situation 
for  a=45°  is  only  slightly  "worse".  But  even  for  a=90°  the  cut-off  dis¬ 
tance  of  1.5s  corresponding  to  43.7%  is  acceptable.  Although  it  may  be 


(4.27) 


>  (4. 27’) 


J 
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deemed  marginal,  it  is  still  somewhat  more  satisfactory  than  the  same 
distance  in  the  case  of  geoid  undulations  (under  similar  circumstances, 

1.5s  for  geoid  undulations  is  associated  with  47.6%).  In  summary,  the 
"worst  situation"  corresponds  to  for  ct=90°  or  to  ni  for  a=0,  with  8=45° 
in  both  cases  (interchanging  a<}>  and  AA_  does  not  matter  here).  In  other 
words,  most  affected  by  a  given  point  mass  are  the  rates  of  change  in  £ 
along  the  E-W  direction  and  in  n  along  the  N-S  direction,  as  modeled  at 
observation  points  located  along  the  NE  (and  NW,  etc.)  direction  from  that 
point  mass. 


I 


I 
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D.  OBSERVATION  MODES  IN  TERMS  OF 
A  DOUBLE  LAYER  OF  POINT  MASSES 


An  analogy  between  the  point-mass  representation  of  the  potential 
and  the  density  layer  representation  has  been  exploited  in  the  past.  If 
also  the  double  density  layer  concept  is  applied  to  the  point-mass  model, 
it  leads  to  the  notion  of  a  dipole.  As  in  the  double  layer  situation,  two 
point  masses  (P.M.)  of  the  same  magnitude  but  opposite  signs  can  be  imagined 
located  along  the  same  normal  to  the  reference  surface.  The  depth,  under¬ 
neath  the  earth's  surface,  of  the  upper  P.M.  is  denoted  by  d  as  in  the 
single  P.M.  approach.  The  depth  of  the  lower,  newly  added  P.M.  is  denoted 
by 

d1  =  d  +  v  , 

where  v  is  the  vertical  separation  between  the  twin  point  masses.  And, 
similar  to  the  approach  with  the  single  P.M.,  the  horizontal  separation 
between  the  neighboring  twin  P.M.  is  denoted  by  s. 

In  analogy  to  the  concept  of  double  layer,  the  minus  mass  (or 
minus  density)  has  no  direct  physical  meaning  in  itself  but  is  a  useful 
mathematical  device.  Although  the  number  of  P.M.  themselves  doubles  with 
respect  to  the  single  P.M.  approach,  the  number  of  parameters  remains  the 
same  due  to  the  condition  that  the  P.M.  along  the  surface  normal  differ 
only  in  sign,  not  in  magnitude.  If  v  and  thus  also  the  depth  of  the  lower 
P.M.  grows  indefinitely,  the  characteristics  of  the  single  P.M.  layer  are 
recovered.  On  the  other  hand,  the  decrease  in  v  leads  to  different 
characteristics  which  may  be  desirable  in  some  respects.  For  example. 


)\ 


if  v  becomes  sufficiently  small  compared  to  s,  such  as  v=  0.112s  used  herein 
(corresponding  to  v  =  50  km  for  s  =  445km),  the  relative  effect  of  a  twin  P.M. 
on  the  modeled  value  decreases  much  more  rapidly  with  distance  than  that  of 
a  single  P.M.  With  even  smaller  v  this  feature  is  further  accentuated. 

The  above  discussion  implies  that  a  given  cut-off  distance  in  the 
double  layer  approach  introduces  fewer  approximations  in  the  adjustment 
model  than  does  the  same  distance  in  the  single  layer  approach  (under  the 
implied  assumption  of  constant  d).  The  cut-off  distance  is  defined  in 
[Blaha,  1980]  as  the  distance  from  an  observation  point  beyond  which  no  P.M. 
parameters  are  considered  in  the  formation  of  the  pertinent  observation 
equation.  If  one  tolerates  about  the  same  level  of  approximation  in  both 
approaches,  one  can  decrease  the  cut-off  distance  in  the  twin  P.M.  model  and 
thereby  alleviate  the  computer  run-time  requirements.  In  the  present  analysis 
the  cut-off  distance  corresponds  to  1.5s;  the  resulting  approximations  in 
the  double  P.M.  approach  will  be  seen  to  be  mostly  inconsequential. 

The  present  development  is  based  in  several  respects  on  [Blaha, 
1980],  abbreviated  here  as  [B] .  The  first  section  below  is  concerned  with 
developing  the  twin  P.M.  model  for  five  observational  (and  prediction)  modes 
as  presented  in  [B] .  The  five  kinds  of  modeled  quantities  are:  geoid  un¬ 
dulations,  gravity  anomalies,  deflections  of  the  vertical  (north  and  east), 
horizontal  gradients  of  gravity  disturbance  (north  and  east),  and  vertical 
gradients  of  qravity  disturbance.  The  subsequent  section  deals  with  the  twin 
P.M.  concept  applied  to  the  new  observational  mode,  the  deflection  rates. 

Since  both  sections  are  intended  to  be  essentially  independent,  minor  over¬ 
laps  occur  in  a  few  instances. 


-62- 


5.1  Adaptation  of  Five  Observational  Modes  to  Twin  Point  Masses 


The  geoid  undulation  at  the  (observation)  point  i  due  to  single 
point  masses  is,  according  to  [B] ,  equation  (4.17): 


N.  = 

i 

(1/G)  l  (l/£..)(kM)j  , 

(5.1) 

where 

Z.  . 

ID 

=  distance  between  the  i-th  (observation) 
point  and  the  j-th  point  mass, 

(kM). 

=  j-th  scaled  point-mass  magnitude. 

=  (R2  +  R2  -  2F..)15  , 

(5.2a) 

Ri 

=  R  -  d  , 

(5.2b) 

F.  . 

ID 

=  R  Rx  cos^i;.  , 

(5.2c) 

cos^  =  sin<j>i  sincjK  +  cos4>i  cos<Jk  cos  (Xi~X_.)  , 

(5. 2d) 

where  R  represents  the  earth's  mean  radius  (6.371  km),  G  is  the  average  value 
of  gravity  on  the  earth'  surface  (980  gals),  and  where  is  the  spherical 
distance  between  the  points  i ,  j  with  the  coordinates  (^.Aj  and  (4k,A^), 
respectively. 

In  the  double  P.M.  approach,  the  second  (deeper)  layer  of  P.M.  is 


located  at  the  sphere  of  radius  R'  from  the  earth's  center,  where 


or 


R‘  =  Rx  +  dRi  , 


dR^  =  -v  . 


(5.3a) 

(5.3b) 


The  other  primed  values  are 


F'  =  R  R'  costi*  , 
ij  1  ij 


V  *  (R2  +  R'2  -  2F1  )**  . 
ij  1  ij 


(5.4a) 

(5.4b) 


The  j-th  scaled  P.M.  corresponding  to  Rj  has  the  magnitude  (kM)!  such  that 


(kM')j  =  -(kM)  . 


(5.5) 


The  value  of  the  geaid  undulation  due  to  the  deeper  point  masses 


N‘  *  (1/G)  I  (l/^.)(kM):  . 
j 

Upon  taking  (5.5)  into  account,  the  value  corresponding  to  double  point  masses 


(at  depths  R1  and  R p  is 


N.  =  N.  +  n:  , 

ill 


i  .e. , 


N.  =  (1/G)  l  (1/8,  -  l/i]  )(kM) 


(5.6) 


i:  13  3 


The  formation  of  observation  equations,  etc.,  could  proceed  in  line  with 
Chapter  4  of  [B]  in  conjunction  with  any  observational  mode.  In  particular, 
a  row  in  the  matrix  A  as  described  in  Section  4.1  of  [B]  could  be  formed 
by  omitting  (kM)_.  and  utilizing  the  indices  j  to  determine  the  ordering  of 
the  elements  in  this  matrix. 
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The  gravity  anomaly  at  i  due  to  the  single  point  masses  is 


Ag.  =  (1/R)  l  (l/£  ,)[(R2-  F  )/£2.  -  2](kM)  ,  (5.7) 

1  .  id  id  13  j 


as  can  be  gathered  from  (4.19)  of  [B] .  In  analogy  to  the  above,  the  gravity 
anomaly  in  the  double  P.M.  approach  is 


Ag.  =  (1/R)  l  {(1/4  )[(R2-F  )/£2  -2]  -  (1/4'  )[(R2-  F1  )  /  4 1 2  -  2]  }  ( kM ) 
1  j  ij  ij  iD  ij  ij  ij 

(5.8) 


The  deflections  of  the  vertical,  £  and  ni»  in  the  single  P.M. 
approach  read  (see  4.24  and  4.25  in  [B]): 

(5.9a) 

E.  =  -(R/G)£  ( l/£3  ) [cos<p  sin<j>  -  sin<|>  cos<p  cos(X  -X  )](kM)  , 

1  Tjij  ij  ij  ij  j 

n  =  (R/G)£  ( 1/ £3 .  )cos<J>  s i n ( A  -X.)(kM).  .  (5.9b) 

i  1  j  iD  j  i  D  D 


The  corresponding  expressions  in  the  double  P.M.  approach  are 

(5.10a) 

E  =  -( 1/G)  l  (R  / £3  -  R'/£'3) [cosd  sin^  -  sin<|>  cos<)>  ccs(x  -X  )] ( kM) 
i  j  1  ij  ij  i  j  i  j  i  j  j 


n.  =  (1/G)  l  (R  /S>3  -R'/£'3)  COS 4)  sin(X  -X  )(kM)  .  (5.10b) 

1  j  1  ij  1  ij  D  1  D  j 


The  horizontal  gradients  of  gravity  disturbance,  Tzx  and  Tzy  , 

i  i 

in  the  single  P.M.  approach  read  (see  4.28  and  4.29  in  [B] ) : 


T  =  -3R.  7  ( 1/ £5  )(R -  R  cosij;  )[cos<J>  sin<j>  -  sincp  cos<t> 
zx.  1  ij  1  ij  i  j  i  j 

«  cos(Xi  -  X ,  )](kM)..  , 


(5.11a) 


(5.11 


T  =  3R.  7  [l/l  .)(R-R.  cosijj.  . )cos(j).  sin(X.  -  X . ) ( kM) .  . 
zy.  1  L  13  1  13  3  1  3  3 

1  j 

In  the  double  P.M.  approach  the  corresponding  expressions  are 

T  =  -3  7  f(R,/Jt5.)(R- R,  cos^.  .)  -  (R'/fi-! 5)(R-  R'  cos^..)] 
zx.  “  L  1  13  1  13  13  1  13  J 

1  j 

*[cosi})i  si n<f>_.  -  sindi  cos^  cos(Xi  -  X  )]  (kM)  j  ,  (5.12, 

V  =  3  ?  I(Ri/Jlij)(R- Ri  cos*ii>  -  r;  ^s*..)] 

»cos<}>.  sin(X.  -  X  )(kM) .  .  (5.121 

3  1  r  d 

The  vertical  qradient  of  gravity  disturbance  at  i  in  the  single  P.M. 
approach  reads  (see  4.31  in  [B]): 

Tzz  =  l  cos^  )2/^  -  l] (kM)  .  (5.13 

i  j  j  j  j  j 

In  the  double  P.M.  approach  the  corresponding  expression  is 

Tzz  =  l  ((1/^.)[3(R- Rx  cos^..)2/*2.-  1]  -  (1/A^)[3(R-RJ  cosij,..)2/* 

-  l]}(kM)  .  (5.14 

Next,  a  decision  regarding  a  cut-off  distance  when  adjusting  a 
double  P.M.  model  will  be  made,  again  along  the  lines  of  [B] ,  Chapter  4. 

In  this  analysis,  one  twin  point  mass  will  be  considered,  (kM)_.  and  (kM)^ 
in  conjunction  with  the  observation  point  i  on  the  earth's  surface  (here  a 
sphere  of  radius  R  =  6,371  km).  The  depth  (d)  of  the  shallower  P.M.  is 
taken,  in  agreement  with  the  previous  chapter,  as  a  0.8-multiple  of  the 


horizontal  separation  (s)  between  the  twin  point  masses.  Thus,  if  the 
horizontal  separation  in  degrees  is  s°=4°,  then 

s  =  445  km  , 
d  =  0.8s  =  350.0  km  , 

R1  =  6,021  km  . 

The  vertical  separation  between  the  twin  P.M.  is  chosen  as 
v  =  50  km  , 

from  which  it  follows  that 
R|  =  5,971  km  . 

The  magnitude  of  (kM) ,  is  chosen  to  be  a  10“6-part  of  the  earth's  kM 
(kM  =  3.986xl014  m3/sec2).  Next,  the  position  of  the  observation  point  i 
is  varied  and  the  corresponding  value  Ni  computed.  This  variation  can  pro¬ 
ceed  by  selected  multiples  of  s,  which  make  it  possible  to  see  at  how  many 
s-intervals  away  from  the  observation  point  the  P.M.  parameters  will  still 
significantly  affect  the  modeled  value  (and  vice  versa)  and  should  thereby 
be  included  in  the  observation  equation  or  in  the  prediction  formulas.  In 
theory  one  would  include  all  of  the  P.M.  parameters  in  every  observation 
equation,  but  this  would  often  result  in  prohibitive  computer  run-time  re¬ 
quirements. 

If  one  single  P.M.  were  considered  (instead  of  a  twin  P.M.),  the 
formulas  qiving  N  ,  Ag.  ,  etc.,  would  be  written  as  they  stand  (see  5.1, 

ii 

5.7,  etc.),  except  that  the  symbol  £  would  be  left  out.  In  general,  any 

j 


-67- 


such  observable  quantity  at  i  can  be  symbolized  by  "0  "  and  the  coefficient 

of  (kM).  can  be  denoted  so  that 

D  ID 

0.  =  C. .(kM).  , 

i  ID  D 

where  C_  is  readily  gathered  from  the  pertinent  formulas.  Considering  the 
twin  point  masses,  one  has 

6.  =  (C  -  C'  }(kM)  . 

ij  ij  j 

If  this  model  is  linearized  --  and  thus  the  analysis  is  concerned  with 
relatively  small  separations  between  the  twin  P.M.  along  the  vertical  --  it 
follows  from  (5.3b)  that 

0.  =  -(9C../3R1)(kM).dR1  =  (3C../3R1)(kM).v  . 

This  expression  will  now  be  examined  in  the  context  of  all  of  the  observa¬ 
tional  modes  listed  at  the  outset. 

For  the  geoid  undulation  fL  ,  the  model  reads 
C. .  =  (l/6)(l/£. .)  . 

ID  i-D 

The  straightforward  differentiation  yields 

3C. ./M  =  (1/G)(1/*. ,)a  .  ,  (5.15 

ID  1  ID  ID 

where 

a  =  -(R  -  R  cos^i  )/i2  ,  (5.16 

ID  1  ID  iD 


and  where  the  angular  separation  between  the  points  i  and  j,  is  taken 
in  terms  of  n-multiples  of  s(in  angular  measure).  The  approximate  results 
for  N  are  presented  as  follows: 

n  =  0  ...  16.6  m  (100%)  n  =  1.5  ...  1.6  m  (10%) 
n  =  0.5  ...  10.1  m  (  61%)  n  =  2  ...  0.7  m  (  4%) 
n=l  ...  4.0  m  (  24%)  n  =  2.5  ...  0.4  m  (  2%) 


At  the  distance  Is  (corresponding  here  to  49)  the  effect  of  the  twin  P.M. 
has  diminished  to  about  24%  of  the  maximum  effect  exercised  directly  under¬ 
neath  the  observation  point  (at  n=0).  The  corresponding  distance  could  be 
adopted  as  a  cut-off  distance  with  regard  to  both  observations  and  predic¬ 
tions.  Clearly,  the  influence  of  the  twin  P.M.  on  the  value  of  decreases 
much  more  rapidly  than  is  the  case  with  the  single  P.M.  approach  summarized 
in  Section  3.1,  where,  for  example,  n=2  still  corresponds  to  an  effect  as 
large  as  37.6%. 


The  model  for  the  gravity  anomaly  reads 


ci3  ■  . 

bi3'  ' 


(5.17) 


which  yields 


ac 


/3R  =  (1/R  )(l/i  )  (a.  ,b. .  +  a..c..  -  r.  ) 
1  1  ij  13  13  13  ij  ij' 


cii  *  =  «v2> . 


r.  ,  =  R  cos \li.  Jl2 .  . 
13  13  13 


(5.18) 

(5.19) 

(5.20) 
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The  approximate  results  for  are 

n  =  0  ...  88.0  mgal  (100%)  n  =  1  ...  0.9  mgal  (  1%) 

n  =  0.5  ...  30.4  mgal  (  34%)  n  =  1.5  ...  -1.9  mgal  (-2%) 


At  the  distance  of  merely  0.5s  (here  2°)  the  percentage  is  already  down  to 
34%  and,  even  more  importantly,  it  decreases  so  rapidly  that  at  n=l  it  is 
1%.  A  cut-off  decision  made  in  favor  of  n=l  would  be  more  than  sufficient. 
By  comparison,  in  the  single  P.M.  approach  n-1  corresponds  to  21.4%. 


Similar  to  [B] ,  the  analysis  of  the  deflections  of  the  vertical 
leads  to  the  same  results  whether  or  n.  is  considered.  The  case  with  n. 

i  i  i 

yields 


(R./G)(l/A?.)s.  .  , 

1  i]  i] 


s . .  =  cos(j> .  sin(X.  -  X . )  , 
ID  D  ID 


(5.21) 


where  the  point  j  is  considered  at  the  equator,  thus  cos4>_.  =  1.  Subsequently, 
one  has 

8C../3R1  =  (R1/G)(l/^.)s..(l/R1  +  3a..)  .  (5.22) 

The  approximate  results  pertaining  to  Hi  (or  ^)  are 

n  =  0  ...  0.0"  (0%)  n  =  1.5  ...  1.3"  (16%) 

n  =  0.5  ...  8.0"  (100%)  n  s  2  ...  0.5"  (  6%) 

n  =  1  ...  3.5"  (  44%)  n  =  2.5  ...  0.2"  (  3%) 


A  cuc-off  distance  at  1.5s  (here  6°)  would  be  more  than  sufficient. 


corresponding  to  about  16%  of  the  maximum  effect  for  n  =  0.5.  By  comparison, 
the  single  P.M.  approach  shows  about  the  same  effect  (more  precisely  17.1%) 
for  n=  3s. 


The  analysis  of  the  horizontal  gradients  of  gravity  disturbance 
produces  the  same  results  for  either  Tzx  or  Tzy  gradient.  The  case  with 
Tzy  yields 


C.  . 
ID 


3R,(l/£: .)u. . 

1  ID  ID 


S.  . 
ID 


» 


u.  .  =  R  -  R,  cos \p.  .  , 

ID  1  riD 


(5.23) 


where  j  is  again  taken  at  the  equator;  thus 


3C.  ,/3R  *  3R.(1/A?.)u.  .s..(l/R.  +  5a.  .- cos<».  ./u..)  .  (5.24) 

ID  1  1'  ID  iD  ID  1  3-D  1J  ID 


o  o 

?he  approximate  results  in  E  (Eotvos,  0.1  mgal/km  or  10  sec  )  pertaining 

to  T  (or  T  )  are 
zy.  zx. 


n  =  0  ...  0.00  E  (  0%)  n  =  1  ...  0.49  E  (14%) 

n  =  0.25  ...  3.45  E  (100%)  n  =  1.5  ...  0.03  E  (  1%) 

n  =  0.5  ...  2.86  E  (  83%)  n  =  2  ...  -0.02  E  (-0.5%) 


A  cut-off  distance  at  Is  (here  4°)  would  be  more  than  sufficient,  correspond¬ 
ing  to  about  14%  of  the  maximum  effect  at  n=  0.25.  By  comparison,  the  single 
P.M.  approach  shows  about  the  same  effect  (more  precisely  17.2%)  for  n=  1.5. 
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The  approximate  results  for  T  are 

i 


n  =  0  ...  7.96  E  (100%)  n  =  1  ...  -0.35  E  (-4%) 

n  =  0.5  ...  1.16  E  (  15%)  n  =  1.5  ...  -0.18  E  (-2%) 


A  cut-off  distance  at  0.5s  (here  2°)  appears  more  than  sufficient,  and  it 
corresponds  to  about  15%  of  the  maximum  effect  at  n=0.  By  comparison,  the 
single  P.M.  approach  shows  a  37.5%  effect  for  the  same  distance. 
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5.2  Development  of  Deflection  Rates  in  Terms  of  Twin  Point  Masses 


The  development  in  this  section  is  based  on  the  results  from 
Section  4.3,  where  the  model  for  deflection  rates  has  been  formulated  in 
terms  of  single  point  masses.  The  notations  of  the  previous  section  are 
adopted  without  changes,  with  primes  indicating  the  quantities  connected 
with  the  deeper  layer  of  the  point  masses,  designating  the  quantities 
resulting  from  the  effect  of  both  layers  of  point  masses,  v  denoting  the 
vertical  separation  of  the  twin  point  masses,  etc.  Again,  the  twin  point 
masses  have  the  same  mass  magnitude  but  the  opposite  sign.  Thus,  from 
(4.18)  and  (4.19)  we  write  the  basic  model  for  the  deflection  rates  as 
follows: 


k.  =  (R./GR)  I  {[(l/fc3.)<p.  .t. .  -  tg<p . )  -  (l/£.,3)(p:.t.  .  -  tg*. )] 

1  1  1]  1]  1]  1  1]  X]  1]  1 


>  (q.  ./cos<J>.  )sina  -  f(  1/e3  . ) (p.  .t2  .  -  cosiJj.  . ) 

^ii  i  L  xi  ii  ii  *11 


"(l/^ij)(Pijtij  ’  C0S^ij)]C0SaXkM)  j  » 


(5.26) 


n  =  (R  /GR  cos4>.)  T  { [( 1/&3  )(p  q2  -  cos<i>  cos<j>  cosAX  .) 
i  1  1  il  ij  ij  i  j  il 


•l  n 

-(1/2.!  )(p!  .q.  .  -  cos*.  cos<p .  cosAX.  .)!  (l/cos<f>.  )sina 

ij  i]  i]  i  j  ii  J  i 


-(p.  ./I2 .  -  p!  ,/2.!3)q.  .t.  .cosa)(kM) .  , 

ii  ii  ii  ii  i]  ii  1 


(5.27) 
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where  the  trigonometric  functions  associated  witn  the  uoeper  point  mass  are 
unchanged  when  passing  to  the  primed  quantities,  including  cosy  ,  t  and 

q .  . ,  but  where 

13 


3  RR*/^  -*2  • 
1  13 


As  in  the  approach  with  the  single  point  masses  in  Section  4.3,  we 
consider  only  the  effect  of  one  P.M.  magnitude  located  at  the  point  j  (where 
<{>  =0)  on  the  model  value  at  the  (observation)  point  i.  Adopting  also  the 
previous  notations  and  approximations,  we  can  derive  --  or  transcribe  from 
(4.20)  and  (4. 20')  --  the  following  relations: 


(5.28) 


i  =  -(R./GR)((p.  ./£3.  -  p:  ./V.  3)A<i>.  .A,\.  .sina  +  [(p.  ./I2  .  -  p!  ./V.  3)Ad>2  . 
1  l  13  13  13  13  13  13  L  13  13  13  13  13 


-  (l/Z*.  -  l/Z!p]cos<x}(kM).  , 


*  -(R^GVfKp.JZ*  -  p'/Z!pA\'  -  (1/Z*  -  1/Z!p]sincc 


(5.28* ) 


+  (Pij/ilij'Plj/)li?)^i3AXi3COSa}(KM)3; 


the  comment  made  in  conjunction  with  (4.20)  and  (4.20* )  can  be  adopted,  word 
for  word,  also  in  the  present  situation  as  depicted  in  the  above  two  equa¬ 
tions.  Accordingly,  only  (5.28)  will  be  examined  with  respect  to  a  and  B. 

With  regard  to  the  vertical  separation,  we  adopt 
v  =  0.112s  , 
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which  corresponds  to  v=50km  for  s  =  445km  (this,  in  turn,  corresponds  to 
4°).  Thus,  recalling  that  in  angular  measure  s  has  been  denoted  as  Aw,  we 
have 

d'  =  d  +  v  =  0.912s  =  0.912  R  Aw,  d'2  =  0.83  R2  Aw2  . 


If  the  value  0.64  in  the  formulas  (4.22)  is  replaced  by  0.83,  it  can  readily 
be  shown  that 


l/£3.  ■  1/£ij  =  (1/r3  Aw3 ) [0 . 285/ ( n2  +  0.64) 5/2]  , 
Pij/£ij'  pij/£ij  *  (1/r3  ^5>  [l.425/(n2+  0.64)7/2]  , 


(5.29) 


provided  0.83  -  0.64  =  0.19  can  be  considered  small  compared  to  n2  +  0.64, 
which  is  quite  acceptable  for  the  values  we  are  most  interested  in  (between 
n=l  and  n=2).  The  relations  in  (5.29)  are  derived  from 

d(l/x)k/2  =  -(k/2)(l/xk/2+1)dx  , 

with  x=n2+  0.64,  dx=  0.19  and  k=  3  or  5. 

If  (5.28)  is  developed  with  the  aid  of  (5.29),  as  well  as  with  the 
notations  introduced  in  (4.21),  it  follows  that 

ki  =  -C[l/(n2+  0.64)7/2] [l.425n2sinB  cosB  sina 

+  (1.425n2cos2S  -  0.285n2  -  0.182)cosa]  ,  (5.30) 


M 


< 


ti 


tl 


where  C  is  the  same  as  in  (4.231).  Since  in  view  of  (4.24),  (4.24')  we 
again  have  a=b,  an  extreme  for  the  function  in  the  second  brackets  above 
does  not  exist,  and  we  proceed  to  examine  a=0,  45°  and  90°,  searching  now 
for  extremes  with  respect  to  8.  The  strategy  used  for  the  single  point 
masses  will  be  pursued  here  as  well. 

Azimuth  g=0.  In  exact  analogy  to  Section  4.3,  a  maximum  occurs 
at  8=0  which  now  coincides  with  an  absolute  maximum  already  for  n>0.65 
(previously  this  held  true  for  n>  1.13).  Parallel  to  (4.25)  and  (4. 25') 
we  have 


a=0,  6=0;  £.  =  -C[(1.140n2  -  0.182)/(nz  +  0.64)7/2]; 


(5.31) 


n 

=  0 

. .  100% 

n  =  1 

..  -19.5%  1 

=0.5  ... 

..  -17.8% 

n  =  1.5  . . . 

..  -6.7% 

>  (5. 31') 

n 

=  0.69  .. 

..  -28.3% 

n  =  2  ... 

..  -2.3% 

Note:  In  addition  to  the  primary  maximum  at  n=0,  a  secondary  maximum 
exists  at  n=0.69.  Beyond  n=l  the  values  in  (5.31')  are  seen  to 
decrease  much  more  rapidly  than  their  counterparts  in  (4.25’ ) . 

Azimuth  g=45°.  A  maximum  occurs  at  3=22.5°  as  in  Section  4.3, 
but  it  now  corresponds  to  an  absolute  maximum  already  for  n>0.56  (as 
opposed  to  n>0.89  found  previously).  In  analogy  to  (4.26)  and  (4. 26')  we 
have 

a=45° ,  3=22.5°;  =  -C[( 1 .435n2  -  0. 182)/(n2 +  0.64)7/2] ,  (5.32) 
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n  =  0 

.  100% 

n  =  1 

_  -25.6% 

n  =  0.5 

.  -30.6% 

n  =  1.5  . 

....  -8.6% 

>  (5.32') 

n  =  0.66 

.  -39.6% 

n  =  2 

_  -3.0% 

Note:  In  addition  to  the  primary  maximum  at  n=0,  a  secondary  maximum  exists 
at  n=0.66.  Beyond  n=l,  the  values  in  (5. 32')  decrease  much  more 
rapidly  than  their  counterparts  in  (4.26'),  similar  to  the  previous 

case. 

Azimuth  a=90°. 

Here  again,  the  (absolute)  maximum  at  B=45( 

3  cor- 

responds  to 

its  single  P.M.  counterpart. 

Similar  to  (4.27)  and  (4.27'), 

the  results 

are 

a=90° , 

6=45° ;  £ . 

=  -C[0.712n2/(n 

2+  0.64)7/2]  ; 

(5.33) 

n  =  0 

.  0 

n  =  1 

....  -47.1%  'J 

n  =  0.5  , 

.  -100% 

n  =  1.5  . 

....  -14.6% 

)  (5. 33') 

n  =  2  .  -4.9%  J 


Note:  The  formula  (5.33)  reveals  only  the  (primary)  maximum  at  n=0.5  (more 
precisely  at  n=0.51) . 

The  case  a=90°  is  again  much  "worse"  than  a=45°  which,  in  turn, 
is  only  slightly  "worse"  than  a=0.  This  is,  accordingly,  the  "worst  situa¬ 
tion"  as  discussed  previously  in  Section  4.3  dealing  with  the  single  point 
masses.  But  even  for  ot=90°  the  cut-off  distance  of  1.5s  is  more  than  suf¬ 
ficient  (it  is  associated  with  14.6%  while  in  the  single  P.M.  approach  we 
had  43.7%).  We  can  again  summarize  this  development  by  stating  that  the 
"worst  situation"  occurs  for  r.  at  a=90°  or  for  n.  at  a=0,  with  3=45°  in 

i  i 

both  cases. 
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6.  CONCLUSIONS 


The  simultaneous  global  adjustment  of  satellite  altimeter  data 
in  terms  of  spherical -harmonic  (S.H.)  potential  coefficients,  tidal  para¬ 
meters  and  state  vector  (s.v.)  components  results  in  a  smooth  approximation 
to  the  oceanic  geoid,  called  the  trend  surface,  and  in  the  residuals  con¬ 
taining  the  suppressed  geoidal  detail.  In  addition  to  the  trend  surface,  the 
adjusted  S.H.  coefficients  can  express  other  features  of  the  smoothed-out 
global  gravity  field.  For  example,  predicted  values  of  gravity  anomalies, 
deflections  of  the  vertical  and  other  functions  of  the  disturbing  potential 
can  be  computed  in  a  desired  grid  from  which  contour  maps  can  be  constructed, 
etc. 

The  satellite  altimeter  adjustments  recently  performed  at  AFGL  are 
based  on  the  short-arc  algorithm,  where  the  arcs'  lengths  have  been  limited 
to  seven  minutes  or  less  and  the  data  consist  of  SEASAT  altimeter  observations. 
Each  satellite  arc  is  described  by  six  s.v.  parameters  considered  independent 
from  arc  to  arc.  The  tidal  parameters  consist  of  a  global  amplitude  factor 
and  a  global  phase  angle  correction  for  each  diurnal  and  semidiurnal  constitu¬ 
ent  included  in  the  adjustment.  Each  long-period  constituent  has  been  at¬ 
tributed  only  the  first  parameter.  The  simultaneous  adjustment  with  an  un¬ 
limited  number  of  arcs  is  feasible  only  by  virtue  of  eliminating  the  s.v. 
parameters  from  the  normal  equations  as  soon  as  the  last  observation  on  a 
given  arc  has  been  processed,  and  of  re-using  the  same  core  space  for  the  next 
arc.  The  short-arc  algorithm  as  well  as  the  weignting  of  altimeter  observa¬ 
tions  and  all  three  groups  of  parameters  have  been  described  in  the  previous 
AFGL  reports. 


The  (weighted)  tidal  parameters  are  shown  in  Chapter  2  to  be  es¬ 
sentially  uncontaminated  by  the  geoidal  errors  or  by  the  systematic  orbital 
errors.  The  insensitivity  of  these  parameters  to  the  geoidal  errors  stems 
from  the  time-dependent  properties  of  the  former  contrasted  to  the  time 
invariabil ity  of  the  latter.  And  the  insensitivity  of  these  parameters  to 
the  systematic  orbital  errors  stems  from  the  global  properties  of  the  former 
contrasted  to  the  local  character  of  the  latter  (the  systematic  errors  change 
from  one  area  to  another).  Such  advantages  do  not  exist  in  a  model  where  the 
adjustment  is  made  in  terms  of  S.H.  tidal  coefficients.  Although  the 
a  priori  values  of  such  coefficients  are  used  in  the  present  model  to  de¬ 
scribe  the  approximate  behavior  of  the  pertinent  diurnal  or  semidiurnal  con¬ 
stituent,  only  two  of  their  very  special  linear  combinations  represented  by 
the  above  two  tidal  parameters  are  subject  to  adjustment.  Clearly,  non- 
adjustable  coefficients  cannot  be  affected  by  orbital  or  other  systematic 
errors. 

Chapter  2  also  illustrates  the  improvements  in  the  trend  surface, 
geoidal  resid-als  and  "observed"  geoid  undulations  due  to  the  inclusion  of 
the  tidal  adjustment  in  the  overall  adjustment  of  satellite  altimetry.  The 
geoidal  residuals  can  be  used  as  observations  in  a  subsequent,  or  second- 
phase,  adjustment  of  a  short-wavelength  oceanic  geoid  in  terms  of  point-mass 
(P.M.)  magnitudes  as  parameters.  And  the  "observed"  geoid  undulations,  ob¬ 
tained  by  superimposing  the  geoidal  residuals  on  the  trend  surface,  can  serve 
in  describing  the  geoidal  detail  in  yet  another  fashion.  As  one  example, 
they  can  be  used  to  produce  a  high-degree  and  order  set  of  S.H.  potential 
coefficients  via  integral  formulas  (not  via  an  adjustment),  which  can  then 


serve  in  predicting  the  desired  geophysical  quantities.  Further  properties 
and  uses  of  the  geoidal  residuals  and  "observed"  geoid  undulations  are  listed 
in  the  latter  part  of  Chapter  2. 

Chapter  3  is  concerned  with  a  reformulation  of  the  second-phase 
adjustment  in  terms  of  P.M.  parameters,  leading  to  a  banded  system  of  normal 
equations.  Such  a  system  can  be  arrived  at  and  resolved  by  addressing  three 
separate  tasks.  The  first  task  is  to  eliminate  a  number  of  point  masses 
from  a  given  observation  equation.  However,  these  point  masses  should  be 
located  sufficiently  far  from  the  pertinent  observation  point  so  that  the 
rigor  of  the  solution  not  be  unduly  compromised.  Due  to  the  bandwidth  con¬ 
siderations,  the  cut-off  distance  should  be  as  small  as  practicable.  This, 
in  turn,  implies  that  the  depth/side  ratio  characterizing  the  P.M.  distribu¬ 
tion  should  be  reasonably  small.  As  a  result,  the  1.6/1  ratio  used  recently 
has  been  lowered  to  a  0.8/1  ratio.  The  second  task  is  to  arrange  the  P.M. 
parameters  in  such  a  way  that  a  banded  structure  may  indeed  materialize. 

The  first  task  alone  would  merely  lead  to  a  great  number  of  zeros  in  the 
normal  equations,  not  necessarily  to  their  banded  structure.  And  the  third 
task  is  to  solve  the  resulting  system  with  a  maximum  efficiency  through  an 
adaptation  of  the  well-known  Choleski  algorithm. 

The  combined  solution  to  the  three  tasks  above  gives  rise  to  a 
"modified  Choleski  algorithm".  Whereas  previously  only  small  areas  could  be 
resolved  in  a  P.M.  approach,  the  oceanic  geoid  over  entire  ocean  basins  can 
now  be  adjusted  in  a  few  overlapping  strips  of  point  masses.  Indeed,  the 
modified  Choleski  algorithm  reduces  both  the  run-time  and  the  core-space 
requirements  several  fold.  For  example,  1,300  point  masses  can  be  deployed 
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with  this  algorithm  as  compared  to  about  200  point  masses  which  would  have 
been  admissible  otherwise. 

In  Chapter  4  a  model  for  the  deflection  rates  along  a  given  route 
is  developed  in  tensor  notations,  and  is  subsequently  specialized  for  the 
S.H.  and  P.M.  parameters.  It  is  pointed  out  that  the  rate  of  change  in  the 
north  deflection  of  the  vertical  (0  when  proceeding  eastward  is  not  the 
same  as  the  rate  of  change  in  the  east  deflection  of  the  vertical  (n)  when 
proceeding  northward.  However,  the  difference  between  these  two  rates  is  ex¬ 
pressed  in  a  simple  relationship  which  can  serve  as  a  useful  verification. 

The  formulas  developed  in  terms  of  the  S.H.  potential  coefficients  are  veri¬ 
fied  in  Chapter  4  and  the  formulas  developed  in  terms  of  the  P.M.  parameters 
are  verified  in  the  Appendix.  Chapter  4  closes  with  an  analysis  of  a  cut¬ 
off  distance  under  varying  circumstances.  It  is  concluded  that  even  in  the 
worst  case  situation  the  approximations  associated  with  a  recommended  cut-off 
distance  are  somewhat  less  serious  than  the  approximations  associated  with 
the  same  distance  in  the  case  of  geoid  undulations. 

Chapter  5  transforms  the  observational  modes  developed  and  described 
in  the  previous  AFGL  reports  from  the  context  of  single  point  masses  to  the 
context  of  twin  point  masses  (the  former  implies  a  single  P.M.  layer  while 
the  latter  implies  a  double  P.M.  layer).  These  modes  are:  geoid  undulations, 
gravity  anomalies,  deflections  of  the  vertical  (north  and  east),  horizontal 
gradients  of  gravity  disturbance  (north  and  east),  and  vertical  gradients  of 
gravity  disturbance.  In  the  latter  part  of  Chapter  5  the  same  transformation 
takes  place  also  for  the  newly  developed  deflection  rates.  The  deflection 
rate  formulas  in  terms  of  twin  P.M.  parameters  can  be  verified  in  close 
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analogy  to  their  single  P.M.  counterparts.  Since  the  same  verification 
takes  place  for  the  upper  and  the  lower  point  mass,  and  since  the  contribu¬ 
tions  of  these  two  point  masses  are  added  algebraically,  the  pertinent  for¬ 
mulas  in  the  twin  P.M.  model  are  automatically  verified.  An  analysis  of  all 
the  observational  modes  considered  confirms  the  plausible  property  that  if  a 
single  point  mass  is  made  into  a  twin  point  mass,  the  effect  exercised  on  a 
modeled  value  decreases  much  more  rapidly  with  distance. 


APPENDIX 


VERIFICATION  OF  THE  POINT-MASS  MODEL  FOR  THE  DEFLECTION  RATES 


The  verification  formula  (4.8)  for  the  deflection  rates,  applied 
at  the  observation  point  i,  reads 


SE  -  nN  =  ( l/GR2cosd>i ) tg<J>i  3T./3A.  , 

i  i 


where  3T./3A.  is  given  in  [Blaha,  1980],  equation  (4.7b),  as 


3T./3A.  =  -RR.cos<}>.  I  [1/1.  .)cos<j>.sinAA.  .(kM) . 

i  i  1  ri  4  13'  r3  13  3 


with  AX.  .  denoting  X. -  X.  as  before.  We  then  have 
13  1  D 

iE  -  \  =  -(R1/GR)tg<f)i  £  (l/£3  )cos<}>.sinAA  (kM) 
i  i  j  j  j  j  j 


(A.  1) 


On  the  other  hand,  the  formula  (4.18)  with  a=90°  and  the  formula 
(4.19)  with  a=0  yield  respectively: 

=  (R  /GRcosq  )  £  (l/£3  )(p  t  -  tg<j>  )q  (kM)  , 

n  =  (R./GRcos<J>. )  l  (l/t3,)p.  ,q.  ,t.  ,(kM) .  , 

'n  1  1  .  13  ri3  13  13  3 

i  3 

and  thus 


CE.-0N_  =  -(R1/GRcos4>.)  I  (l/2-3j)tg4>iqij(kM)_. 


(A. 2) 
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But  due  to  the  notation  introduced  following  (4.19),  namely 


q.  .  =  cosdj.cosA.sinAA.  .  , 
3-D  1  D  ID 


equation  (A. 2)  becomes  (A.l)  and  the  verification  is  terminated. 
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